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The present version is a further development of the program described in 
Ref. 7 on page i. In addition to the contributors already named on page i, F. H. 
Whitlock has contributed materially to the programming and the checkout of the 
additional modifications. The authors further wish to thank Mrs. A. Michalowski 
and Mrs. F. Magee for their assistance in the preparation of this manual. 

It might be useful at this point to draw attention to some major changes 
effected in update 12C . They may be summarized as follows: 

1 . The Kepler problem has been completely reformulated in terms 
of universal variables and has been programmed in hardware 
double precision. This formulation is described in Appendix D. 

2. The position and velocity vectors are carried in double precision, 
even though the perturbations are still computed in single pre- 
cision only. 

3 . The perturbation acceleration computation has been changed to 
allow for the inclusion of additional planets as well as for addi- 
tional harmonics of the earth and moon gravitational fields . 

4. The planetary ephemerides are read from tapes derived either 
from cards supplied by the Naval Observatory or from the JPL 
ephemeris tapes . The coordinate system in each case is chosen 
as the mean equator and equinox of the middle of each two-year 
file. 

5 . The integration routine has been changed to permit integration 
in the negative time direction and the use of the universal anomaly 

as independent variable . The computation of the typical term 
3 3 

R/r - R 0 /r o in the planetary perturbations has been reformu- 
lated to avoid both the loss in accuracy and the binomial series 
expansion, as described in Appendix E . 

6. A new trajectory search routine has been incorporated, using 
the techniques described in Reference 6. For this routine, the 
number of dependent and independent variables need not be equal. 

It is described in Section IV-H and Section vm-C-1. 

7. The integration and print interval control has been generalized 
and combined with the termination logic. It is described in 
Section VIII -G. 
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8. A new shadow routine has been added which is described in 
Section VIII -X and Appendix N. 

9 . The addresses on the modification cards are to be given in 
octal form, thus making laborious conversions unnecessary. 

10, All cards read into the program (modifications as well as 

regular input) are printed, thus a record of all the parameters 
introduced into a run is kept. 
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I. 


INTRODUCTION 


This report describes a general purpose Interplanetary Trajectory 
Encke Method (ITEM) Program, programmed for the IBM 7090 and IBM 
7094. The method employed is designed to give the maximum available 
accuracy without incurring prohibitive penalties in machine time. On the 
basis of research described in Reference 4, the Encke method was selected 
as best satisfying these requirements . However, the classical Encke 
method was modified to eliminate some of its objectionable features. This 
modified Encke method is described in Appendix A. 

The perturbations included in this program are the gravitational 
attractions of earth, moon, sun, Jupiter, Venus and Mars considered as 
point masses. Additionally, the effects of the second, third and fourth 
zonal harmonics and the tesseral harmonics of the earth and moon gravi- 
tational fields, as well as the aerodynamic drag, small corrective thrusts 
and radiation pressure including the shadow effect of the earth, are con- 
sidered. The input may be prepared in any one of several common sys- 
tems and a great variety of output options is available . 
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II. NOTATION 


Upper case - vectors; Hats - unit vectors; Lower case - magnitudes 


Description 

Symbol 

Units 

Cartesian coordinates of vehicle with respect 
to reference body 

x y z 

km 

Velocity components in Cartesian coordinates 

x y z 

km/sec 

Time 

t 

hrs. 

Longitude measured from Greenwich, + East 
(used in Section IV and Appendix H) 

e 

degrees 

Longitude of vernal equinox 

6 o , 

degrees 

Speed 

V 

km/ sec 

Geodetic altitude* 

h 

km 

Geodetic latitude 

4> 

degrees 

Geodetic flight path angle 

7 

degrees 

Geodetic flight path azimuth 

A 

degrees 

Acceleration parameter (defined in Appendix E) 

u 


Right ascension 

RA 

degrees 

Astronomical units 

AU 


Earth radii 

ER 


Earth mass 

m 

e 



*Note: The following 3 symbols with primes denote the corresponding geocentric quantities. 
Geocentric in this report refers to a spherical earth, L e., e 2 = 0 . In this case 
= 8 - declination. 



Description 

Vehicle position vector 

Distance to vehicle 

Perturbation displacement vector 

Perturbation displacement vector components 

Perturbation acceleration 

Coordinate functions and their time derivatives 
respectively 

Mass parameter 

Semi-major axis 

Earth's eccentricity as used in Appendices 
H , I, L, S 

Mean motion 

Unit vectors for classical two-body orbit solution 
Eccentricity anomaly as used in Appendix D 
Elevation angle as used in Appendix I 

*o * R o 

Incremental eccentric anomaly as used in 
Appendix D 

Functions of 6 defined by equations (E.2) 

Inclination of orbital plane 

Right ascension of the ascending node 


Symbol 

R 

r 

AR 

£> t ?> l 

F 

f. g. f. g 


e 

n 

XV XV 

p, Q 

E 

E 

d o 

e = e - e 0 

f f f f f 

0 ’ 1 ’ 2 ’ 3 ’ 4 


fl 
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Description 


Symbol 


Argument of perigee &> 

Parameters which account for polar oblateness of 

the earth, defined in Appendix H c , s 

Right ascension of the station meridian RA s 

Range measured from observation station p 

Direction cosines measured in a topocentric 

coordinate system A., p, v 

Declination S 


SUBSCRIPTS 

Vehicle 

ith perturbing body 

Quantity obtained from Keplerian solution of 
two-body problem 

Reference body as used in Appendix B 
Station 

r a - R b 

Value at rectification time 

Corresponds to x, y, z components respectively 
Value at perigee 


v 

i 


k 

c 


s 



o 


n = 1, 2, 3 


P 
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HI . GE NERAL PROCEDURE FOR USING PROGRAMS 

Initial conditions, terminal conditions and print frequency, as well as 
other parameters controlling the flow of the program, are read as input. The 
computation of the trajectory then proceeds until one of the terminal conditions 
(e.g., maximum time) has been reached or an error is encountered. At this 
time the program prints the reason for its termination and selected storage 
locations, then proceeds to the next case. When an end of file is encountered 
on the input tape, control is transferred to the monitor, 

A modification card which causes a core dump to be printed is available . 
An option is also included which permits the program to iterate on initial con- 
ditions for trajectories in order to satisfy a set of desired conditions at the 
moon, earth, or a designated planet. 
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IV. 


INITIAL CONDITIONS 


The initial conditions necessary for the specification of a trajectory 

are: 

1. Initial position of the vehicle relative to the reference body. 

2 . Initial velocity of the vehicle relative to the reference body. 

3 . Initial time of launch referenced to a base time . 

For specification of the initial conditions, the reference systems and 
units shown below may be used. 

A. Cartesian Coordinates 

The coordinate system is defined as follows: 

1 . The origin is at the center of the reference body. 

.2. The x-axis is in the direction of the mean equinox of December 31.0 
the year specified as year of launch. 

3. The xy-plane is the mean equatorial plane of the earth. 

Initial position is given by the x, y, z coordinates of the vehicle. 
Initial velocity is given by the x, y, z components of the vehicle . 
Initial time of launch from base time 1 (t) is also given. If the 
program is used in its standard form, the units ^ to be used for 
the above are: 

x , y, z - km 

x , y, z - km/sec 
3 

t - days, hours, minutes, seconds from base time. 

The year of launch must also be given. 


1 The base time is 0. Oh U. T. December 31 of the year previous to the year of 
launch. 

2 Scale factors are used to convert the input units to the units used internally 
(ER or AU and hrs) . Any other set of units may be used by changing these 
scale factors with a modification card as described in Section VTH. 

3 The number of days must be an integer. Hours are expressed as floating point 
decimals; minutes and seconds may be included in the number of hours as a 
fractional part. Zeros must then be loaded in place of minutes and seconds. 
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B. Geodetic Polar Coordinates 

Initial position of the vehicle is given by: 

1. Geodetic latitude (<p) 

4 

2. Longitude, measured from Greenwich (0) 

3. Geodetic altitude (h) 

4 

4. Longitude of vernal equinox at initial time (6^). This quantity 

may be computed by the program or may be loaded. (See IX-A Note 15) 


Initial velocity of the vehicle is given by: 

1. Speed (v) with respect to the center of the earth. 

2 . Flight path azimuth (A) measured clockwise from north in a plane 
normal to the geodetic altitude. 

3 . Flight path angle ( y ) measured from a plane normal to the geodetic 
altitude . 


Initial time of launch from base time (t) must also be given. * 
The following units must be used with the above: 

1. <p, 0, and 0^ - degrees; h - km. 

2. A and y - degrees; v - km/sec . 

3. t - days, hours, minutes and seconds. 


C. Geocentric Polar Coordinates 


Ordinarily an input given in polar coordinates will be interpreted as 
described in paragraph B (preceding). However, if the eccentricity of the 
earth is set to zero, the program will interpret latitude as declination, 


4 If the right ascension (RA) at initial time is known, it may be used in 
place of longitude (0 ). The longitude of the vernal equinox (0 q ) is then 
set to zero . 
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height as distance above a spherical earth of equatorial radius, and flight path 
angle and azimuth with reference to a plane normal to the radius vector. If 
this interpretation of the input is desired, the eccentricity of the earth ( e^ ) 
is set tp zero with a modification card described in Section VTH-A. This ec- 
centricity is used only in this part of the program. Therefore, other portions 
of the program will not be disturbed. An option is also included which inter- 
prets latitude and height as geodetic, and azimuth and flight path angle as geo- 
centric. (See VUT-B-1) 

D . Polar Coordinates in Moon Reference 

If polar coordinate input is used in moon reference, the input quantities 
are interpreted in a righthanded coordinate system defined as follows: 

1. The xy-plane is the moon’s orbital plane . 

2. The x-axis points toward the earth. 

3. The z-axis is parallel to the angular momentum vector of the 

moon about the earth. In this system, 

latitude tp is 

the angle between radius vector and xy-plane, 
positive for positive z . 

longitude 0 is 

the angle between the projection of the radius vector 
on the xy-plane and the x-axis, positive for 
positive y. 

0 irrelevant 
o 

altitude h is 

the distance from the surface of a spherical moon. 

speed v is 

the speed with respect to the center of the moon in 
a non-rotating coordinate system. 

Flight path azimuth and angle are defined in a plane normal to the 
radius vector. 
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E. 


Translunar Plane 


The translunar plane input (see Figure 1) uses the initial conditions 

5 

(height, h, speed, v, argument of radius, and flight path angle, y) 

with respect to the translunar plane . The inclination i and the lead 
6 TX* 

angle cp serve to position the translunar plane with respect to the moon's 
orbital plane. The moon's orbital plane is related to the standard coordinate 
system by means of the ephemeris as described in Appendix W . 

The translunar plane (see Figure 1) input consists of the following: 


Load - 1 
Load 
Load + 1 
Load + 2 
Load + 3 
Load + 4 
Load + 5 


Lead angle (<p) - degrees 
Case No, 

Argument of radius ( "S') - degrees 

Altitude - km 

Speed - km/sec 

Inclination of plane - degrees 

Flight path angle - degrees 


The remainder of the load is the normal polar load. 


5 The angle between the initial position vector and the ascending node with 
respect to the moon’s orbital plane. 

6 The lead angle is the angle between the position vector of the moon at 
launch and the descending node of the translunar plane with respect to 
the moon's orbital plane . 
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Program symbol : 

4> — Lunar lead angle 
i = Inclination of translunar plane 
W = In-plane angle of radius = argument 
of radius 

y = Flight path angle, angle above local 
horizontal (geo or lunicentric) 
h = Injection altitude 


Figure 1. Translunar Plane 



F. 


Ecliptic Coordinates 


Input in ecliptic coordinates is available . The coordinate system used 
in this case is identical to the coordinate system described in Section IV-A, 
with the ecliptic plane replacing the equatorial plane as xy-plane. 


G. Osculating Element Input 

One may input osculating elements as follows: 


LOAD 1 
LOAD - 1 


LOAD 
LOAD + 1 
LOAD + 2 
LOAD +3 
LOAD +4 
LOAD + 5 
LOAD + 6 


4 

Indicator for type of input 

1 = Time of perigee 

2 = Delta - M 

3 = E 

4 — True anomaly 
Case No. 

Argument of perigee 
Longitude of ascending node 
Inclination 

Semi -major axis (in earth radii) 

Eccentricity 

Time of perigee, mean anomaly, eccentric anomaly, or 
true anomaly (depends on whether LOAD -1 has a 1, 2, 

3 or 4, respectively) . 


The program converts the above to Cartesian coordinates and then 
continues normally. 
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H. Trajectory Search 

The trajectory search can determine the initial conditions (in polar 
form) necessary to achieve any set of constraints chosen from among impact 
parameters and osculating elements at a target or a return reference body. 
Its activation is described in VUI-C-1. 

The partial derivatives are computed by the secant method, thus, in 
general, a large number of trajectories must be computed. In the initial 
stages of the search it is recommended that the patched conic mode of com- 
putation be used. Further, throughout the iteration, the generalized anomaly 
should be used as independent variable . The method used in this search is 
described in Reference 6 . 


I. Comments 


1. The program computes in Cartesian coordinates. 

The units used internally in the computation are: 

a. position: Earth radii (ER) - Astronomical units (AU) 

b . velocity: ER/hr - AU/hr 

c. time: hours 

(Earth radii units are used in the earth or moon reference . 
Astronomical units are used in the sun or planet reference .) 

2 . Cartesian or polar coordinates can be used when launching 
from any body, providing ecliptic polar coordinates are 
used in sun or planet reference . 

3 . Equations for converting the initial conditions from polar 
to Cartesian coordinates are shown in Appendix H. 
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V. TERMINATING CONDITIONS 


The set of conditions which will terminate a trajectory may be sum- 
marized as: 

1. Maximum time of flight - hrs. 

2. Maximum distance from any possible reference body considered 
in the solution. 

3 . Minimum distance from any possible reference body considered 
in the solution. 

Any of these conditions will terminate a trajectory. Loading a large 
number into any of the maxima and a zer o into any of the minima will make 
the corresponding conditions inoperative. A proper choice of these numbers 
will permit complete computation of the desired trajectory, avoid extensive 
unnecessary computation and guard against faulty input. 
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V. MOON'S ORBITAL PLANE INPUT AND OUTPUT 


A polar coordinate system is available for input and output which uses 
as its reference plane the moon's orbital plane and the vector from moon 
to earth as unit vector . Polar coordinates in this system are defined anal- 
ogous to geocentric polar coordinates. The cartesian coordinates in this 
system are computed by equations (H.3) with 



and 


e 


0 


0 


Here r B is the radius of the body of departure (earth or moon) . 

These coordinates are then transformed to equatorial coordinates by a 
matrix C computed as follows: 

~ _ R EM 

r EM 

£ - R EM X R EM 

R EM X R EM /v 1 


j .= k 


X 1 
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The transformation matrix C is then given by 



(V .2) 


and 


^ ^ ^MOP 


(V.3) 


® ^ ®MOP 


The matrix C is unitary, and C ' 1 = C*, permitting easy inversion of equa- 
tions (V.2) . 
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yi. 


PERMISSIBLE PERTURBATIONS 


The trajectory computation consists of two parts, the exact solution 
to the two -body problem and integrated additions to this solution for the effect 
of perturbations . The successful control of round-off errors in the modified 
Encke method depends on preventing the accumulated round-off error in the 
integrated perturbation displacement from affecting the computed position. 
This is achieved by always keeping the perturbation displacement small by 
rectifying whenever the perturbation exceeds specified limits. The constants 
mentioned below are used in determining the allowable limits as ratios of the 
perturbation position and velocity to the two-body position and velocity, re- 
spectively. 

This ratio is shown for the position vector in Figure 2. 



Figure 2. Encke Method 
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The recommended values for these ratios are as follows: 


*The parameter u 





is defined in Appendix E. 
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VH. RADAR INFORMATION PROGRAM 


The program may be used to simulate radar data if desired. A maxi- 
mum of 30 stations can be handled at one time. The following information 
is required for each station considered: 

1 . Station name - for identification purposes 

2. Position of radar station 

a. Longitude (6) of the station from Greenwich - 
degrees, minutes, seconds - positive eastward* 

b. Geodetic latitude (<p) of the station - 
degrees, minutes, seconds - positive north* 

c. Altitude (h) of station above sea level - feet 

The simulated radar information consists of azimuth, elevation, topo- 
centric right ascension and declination, slant range and range rate. It is 
printed at every normal print time for which the elevation angle is positive. 
Refraction is not considered. 

This part is coded as a subroutine and may be called at other times, 
if desired. 

In addition, the radar program will now print topocentric right ascen- 
sion and declination as well as Minitrack direction cosines . A mode is avail- 
able which differences the range and Minitrack direction cosines for two 
consecutive runs and prints a summary of these differences. The effect of 
various parameters on these measurements can thus be conveniently studied. 


* Alternatively, these quantities may be given in degrees and decimals . 

0's must be loaded into the positions reserved for minutes and seconds. 

The fractional parts will not appear in the printout reproducing the station 
coordinates. They will, however, be included in the computation. 
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VHI. MODIFICATION CARDS 


Normally, the program proceeds in the fashion described in Section IX. 
Modification cards may be used to operate the program in non-standard fashion 
to provide for special requirements . Any modification card included in a case 
will be operative for all succeeding cases in the. input stack, unless it is re- 
voked . Modification cards are listed in this report with their octal and symbolic 
locations. The octal locations refer to the "update 12C" version of the program, 
and some may be different for later update versions . If such a later version is 
used, the octal locations and the addresses should be checked in the symbol 
table at the end of the assembly listing. The printout contains the update num- 
ber used, following the title print. The octal location is to be punched in col- 
umns 2-6, columns 8-10 contain one of the following operation codes 

DEC, OCT, or BCD 

A number to be loaded must begin in column 12 for DEC 1 and OCT cards . 

The first blank encountered to the right of column 12 terminates reading. A 

2 

BCD card must have a word count in column 12, blanks will not terminate 
reading of BCD cards . 

Some modification cards which are frequently used are described below: 

A. Geocentric Polar Coordinates (i.e., spherical earth) 

Modification card 

30111 IECC DEC 0* 


1 A DEC format will read in fixed or floating point numbers. A floating 
point conversion is performed if the number includes a decimal point 
or the symbol E , 

2 The number of words is the number of symbols (including blanks) divided 
by 6. A word count of 10 is assumed if column 12 is blank. 

* Fixed point numbers 
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To revoke modification card: 


30111 IECC DEC .006693422** 


B. Cartesian Scale Factor - Input and Output 


30112 

REKM 

DEC ' 

(1) 

30025 

MDKM 

DEC .... 

(2) 

30026 

TSCL 

DEC 

(3) 


where 

(1) Number of distance units/ER ** 

(2) 23454.87 x (1) = Number of distance units/AU ** 

(3) Number of time units/hour (used only if time ** 

scale is modified) 

Example: Initial conditions to be in ER and ER/hr. Then card reads: 


30112 RE KM DEC 1. 

30025 MDKM DEC 23454.87, 1. 

To revoke modification card: 

30112 REKM DEC 6378.165 

30025 MDKM DEC 14. 9599 E 7, 3600 


** (1) 

** (2), (3), (4) 

** 

** 


(4) Several numbers appearing on one card, separated by commas, 
are loaded into consecutive locations. 


C. To Initiate Radar 

Modification cards (2 cards) 

12654 RADSTA OCT 007400423630 TSX Z$BGN1, 4 

The above card activates the loading of station coordinates 


** Floating point numbers 
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14771 RADENT 


OCT 007400424002 


TSX ZRADAR, 4 


The above card activates the computation of radar information. To revoke 
modification cards (2 cards): 

12654 RADSTA OCT 076100000000 NOP (1) 

The above card inactivates the loading 

(1) This card must be included with the first repeat case if identical 
station locations are used. 

14771 RADENT OCT 076100000000 NOP 

The above card inactivates the computation. 


D. Table Print - Prints astronomical tables in core at every rectification. 
Modification card: 

13447 TBPRNT OCT 007400425153 TSX TAPRN, 4 

To revoke modification card: 

13447 TBPRNT OCT 076100000000 NOP 

E . Table Print Exits 


Normally the program continues after Table Print 
Modification cards: 

If desired to stop: 


25161 TARTN 

DEC 0 

HTR 

To do next case: 



25161 TARTN 

OCT 002000020310 

TRA SWT 5 

To revoke modification card: 


25161 TARTN 

OCT 002000400001 

TRA 1,4 
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F. 


End Card 


Normally the program returns control to the Monitor after all the 
cases are completed. If a program stop is desired, the following card is 
added after the input. 

37776 SWEF OCT 076100000000 NOP 


G. Integration and Print Intervals are controlled by distances from 
reference body. Seven different intervals are possible. 

In earth reference the following is normal: 

30176 EKVEC DEC 1. , 4. , 125. , 0, 0, 0, 0, 0 Radius vectors** 

30206 EIVEC DEC . 03125, . 5, 0, 0, 0, 0, 0, 0 Integration intervals** 

30216 EPVEC DEC 1. , 6. , 0, 0, 0, 0, 0, 0 Print intervals** 

Between 1 and 4 earth radii, the integration interval will be .03125 

hours and the print interval will be 1. hour. Between 4 and 125 earth radii, 
the integration interval will be . 5 of an hour and the print interval will be 6 
hours . 

In moon reference the following is normal: 


30226 

MRVEC 

DEC .27, 10., 0,0, 0,0, 0,0 

** 

30236 

MIVEC 

DEC . 0625, 0, 0, 0, 0, 0, 0, 0 

** 

30246 

MPVEC 

DEC 1. , 0, 0, 0, 0, 0, 0, 0 

** 


In sun reference the following is normal: 


30256 

SRVEC 

DEC 

.1, 10., 0, 0, 0, 0, 0, 0 

** 

30266 

SIVEC 

DEC 

24, , 0, 0, 0, 0, 0, 0, 0 

** 

30276 

SPVEC 

DEC 

48., 0, 0, 0, 0, 0, 0, 0 

** 


In a planetary reference the following vectors should be read in: 


30306 PRVEC 

DEC (2 to 8 values) 

** 

30316 PIVEC 

DEC (1 to 7 values) 

** 

30326 PPVEC 

DEC (1 to 7 values) 

** 

** Floating point numbers 
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ERVEC and MRVEC should be in earth radii units . SRVEC and PRVEC 
should be in astronomical units . EIVEC, EPVEC, MIVEC, MPVEC, SIVEC, 
SPVEC, PIVEC and PPVEC should be in hours or fractions of ah hour. 


H. Inclusion or Exclusion of Perturbations 

1. Ordinarily included are the gravitational field of the earth (2nd, 
3rd and 4th zonal harmonics) and the gravitational attractions of 
sun, moon, Mars, Venus and Jupiter. To exclude any or all of 
these perturbations, put 0 into the corresponding locations as 


follows; 



30067 

MSI 

DEC 

Sun** 

30066 

MM1 

DEC . . . . 

Moon** 

30072 

MJU1 

DEC 

Jupiter** 

30071 

MMR1 

DEC .... 

Mars** 

30070 

MVE1 

DEC .... 

Venus** 

30101 

J20 

DEC .... 

coefficient of 2nd 
zonal harmonic** 

30102 

J30 

DEC 

coefficient of 3rd 
zonal harmonic** 

30103 

J40 

DEC .... 

coefficient of 4th 
zonal harmonic** 


2 . The perturbation due to the moon's oblateness is computed in moon 
reference. If it is to be omitted, put 0 into 

30142 CONSC DEC 0 * 

3 . The following numbers are to be loaded to restore standard 
operation of the program: 


30067 

MSI 

DEC 

332951.3 

** 

30066 

MM1 

DEC 

. 01229483 

** 

30072 

MJU1 

DEC 

317.887 

** 

30071 

MMR1 

DEC 

, 1078210 

** 

30070 

MVE1 

DEC 

.8147689 

** 

30101 

J20 

DEC 

. 0010823 

** 


* Fixed point numbers 
** Floating point numbers 
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30102 

J30 

DEC 

-2.3 E -6 

** 

30103 

J40 

DEC 

-1.8 E -6 

** 

30142 

CONSC 

DEC 

. 36366998 E -2 

** 


The following perturbations are not ordinarily included: 

4. Radiation Pressure may be included by loading a coefficient into 

30105 RACOE DEC ** 

The number to be loaded is: 

KC A 
r 

m 

2 

C is the radiation pressure in dynes/ cm at a distance of 1 AU 
from the sun. 

(c = 4. 6x 10 5 j (Estimated value) 

cm 

A • 2 

A area in cm 
m mass in grams 

K scaler 360C?(23455 .) 2 /6378. 165 x 10 5 = .11178 xlO 8 
sec to hrs, ER to AU, cm to ER 


The radiation pressure will only be active if sunlight impinges on 
the vehicle . For correct results, the radiation pressure should 
therefore be run only in conjunction with the optional shadow com- 
putation as described in Appendix O . 

If, however, the expected trajectory may be safely assumed to be 
entirely out of the earth's shadow, shadow testing may be avoided, 
with a consequent saving in machine time. In this case, the follow- 
ing modification card must also be included: 

33047 SHADIN DEC 1. ** 


5. If Inclusion of the Aerodynamic Drag is Desired , the drag parameter 
1/2 A/m must be loaded in location COEFL by means of the 

following card: 


** Floating point numbers 
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30104 COEFL 


DEC 


** 

The units of A/m are the area in cnT and the mass in grams. 

A layered atmosphere rotating with the earth is assumed. The 
density is obtained by a linear interpolation of density-altitude 
tables . 

6. To include tesseral harmonics, load 


30153 

TESSTR 

DEC 1 

* 

Normal tesseral coefficients are as follows: 


30114 

J22 

DEC -1.9 E -6 

** 

J 22 

30115 

L22 

DEC -21. 

A 22 

30116 

J31 

DEC -1.51E-6 

-r 

'31 

30117 

1.31 

DEC 0, 

\ * ** 
A 31 

30120 

J33 

DEC - . 149 E -6 

-r 

J 33 

30121 

L33 

DEC 22.8 

\ ** 
A 33 

Moon oblateness coefficients are normally: 


30106 

JM20 

DEC 15. 378959 E -6 

T 

J 20 

30107 

JM30 

DEC 0. 

T ** 

J 30 

30110 

JM40 

DEC 0. 

T 

J 40 

30127 

JM31 

DEC 0. 

T ** 
'31 

30130 

LM31 

DEC 0. 

\ ** 
A 31 

30131 

JM33 

DEC 0. 

T 

J 33 

30132 

LM33 

DEC 0. 

\ ** 
A 33 

30125 

JM22 

DEC 0. 

T ** 
J 22 

30126 

LM22 

DEC 0. 

A 22 


* Fixed point numbers 

** Floating point numbers 



7. If any of these perturbations are to be omitted in repeat cases, 
zero must be loaded in the corresponding locations . 

In order to provide a permanent record of the perturbations in- 
cluded in each particular trajectory, the numbers appearing in 
the locations described in Section VIH-H are printed at the be- 
ginning of each trajectory. 


I. Atmospheric Tables for the Drag Computation are stored in core. 

They are contained in the COSPAR Report accepted for the April 1961 
Meeting in Florence, Italy, "CIRA 1961." If it is desired to change 
this atmosphere, the following procedure has to be followed: 

Modification cards 

33747 NTAR DEC .... the number to be en- 

tered is N-l, where 
N is the number of 
entries in the density 
table . * 

33750 RTBL DEC , the values of r at 

which the density is 
given in ascending 
order, a maximum 
of 54.** 

3 

The units for the air density are grams/cm and the height is given 
in ER from the center of the earth. 

34037 RHOT DEC the values of the air 

3 

density in gm/cm in 
respective order cor- 
responding to the r 
table.** 

If other units are used for the density table, the drag parameter de- 
scribed in Part H of this section has to be read in the same units and 
the constant (-6378. 165 E 5)** normally in 

27763 DRSC DEC .... has to be changed 

accordingly. ** 

The negative sign directs the drag force opposite to the velocity. This 
constant converts the drag from the units used for A, m, and p to ER/hr . 


* Fixed point numbers 
** Floating point numbers 
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J. • Angle from Ascending Node, to Vehicle 
To print angle add the following card: 


26017 EAFANS 

OCT 

002000034701 

TRA ASCSAT 

To revoke: 




26017 EAFANS 

OCT 

076100000000 

NOP 


K . The Program Provides a Special Printout described in the output 
(Section X) near the earth, moon, sun and planets . This printout 
occurs at every integration step and is useful to observe the behavior 
of these relevant quantities during ascent and re-entry. This feature 
is triggered by the following modification cards: 


30336 

SERE 

DEC .... 

earth reference (ER)** 

30337 

SMRE 

DEC 

moon reference (ER)** 

30340 

SSRE 

DEC 

sun reference (AU)** 

30341 

STRE 

DEC * • • 

planet reference (AU) ** 


The numbers used above are the radial distance within which the 
special print is to be effective. The units are Earth Radii for the 
earth and moon references and Astronomical Units for sun and 
planet references . A zero in any of the above locations suppresses 
this feature. 

L. The Rectification Print may be eliminated with the following cards: 

13607 RCTPR+4 OCT 076100000000,076100000000 (1) 

13627 RECTPS OCT 076100000000, 076100000000 (2) 

In-order to obtain the rectification prints in repeat cases, after they 
have been eliminated in one case, the original contents, as given in 
the assembly listing of these locations have to be restored. 

(1) Eliminates ratio print 

(2) Eliminates normal rectification print 


** Floating point numbers 



In order to get the normal print at t = 0 (which is also a rectification), 
the program tests the time. If zero, it prints normal output, otherwise 
the skip print routine is effective . 


M. To Compute and Print Latitude, Longitude and Velocity at Perigee 
add the following card: 

25741 PGCOOV OCT 076100000000 NOP 


N. Print at Ascending Node, Apogee and Perigee 

To obtain the normal print at the crossing of the ascending node, 
insert a fixed point number in: 


30342 

CCNT 

DEC 

Ascending node print* 

44214 

CANT 

DEC .... 

Apogee print* 

44215 

CPNT 

DEC . . . . 

Perigee print* 


The number determines how often the print operates, i.e., 2 means 
print every second crossing of the ascending node. The print occurs 
at end of the first integration interval following the crossing.. 


O. Ordinarily the Printout Includes the Following : 

Time, position with respect to reference body. 

Velocity with respect to reference body. 

More printout may be called for by putting a number (fixed point) into 


33710 KMVEC 
1 in KMVEC +2 adds 
1 in KMVEC +3 adds 
1 in KMVEC +4 adds 
1 in KMVEC +5 adds 
1 in KMVEC +6 adds 
1 in KMVEC +7 adds 
1 in KMVEC +8 adds 


DE C 1 , 1, 

vehicle position with 
vehicle position with 
moon's position with 
vehicle position with 
vehicle position with 
vehicle position with 
vehicle position with 


respect to earth 
respect to moon 
respect to earth 
respect to sun 
respect to Venus 
respect to Mars 
respect to Jupiter 


9 • • • 9 


* Fixed point numbers 



1 in KMVEC+9 adds perturbation displacement 
1 in KMVEC+10 adds perturbation velocity 
1 in KMVEC+11 adds perturbation acceleration 

In addition, printout will always contain 

14517 LLDEC+2 RA and declination 

14715 EAZCAL-5 Coordinates of subsatellite point 

14760 1CONV+5 Moon longitude and latitude in moon reference 

14765 PREXT Osculating elements 

These prints may be eliminated with the following cards: 

14517 LLDEC+2 OCT 002000014521 TRA *+2 

14715 EAZCAL-5 OCT 002000014717 TRA *+2 

14760 1CONV+5 OCT 002000014762 TRA *+2 

To restore: 

14517 LLDEC+2 OCT 007400421821 TSX LLACP, 4 

14715 EAZCAL-5 OCT 007400421300 TSX PPRINT, 4 

14760 1CONV+5 OCT 007400421342 TSX LLMNP, 4 


P . Osculating Elements Print is normally given at each print time . 

To revoke: 

14765 PREXT OCT 002000014770 TRA MOOS CO 

To restore: 


14765 PREXT OCT 007400425163 TSX ELCO, 4 


The Normal Output Refers the Osculating Elements to the Equatorial 
Plane with the x-axis along the mean equinox of date. 


The modification described below will trigger, in addition, the printing 
of the inclination, ascending node and argument of perigee with respect 
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to the moon's orbital plane and the x-axis along the moon-earth 
vector, either at the running time t (see Figure 3) or at an arbit- 
rary fixed time t^ . In addition, this modification will print the 

instantaneous osculating elements of the moon's orbit with respect 
to the equatorial plane, 

26015 ELPT OCT 002000034153 TRA INSMOO 

To use orbital plane at time t^: 

26015 ELPT OCT 002000034127 TRA TOOSC+1 

34126 TOOSC DEC time (1)** 

(1) If no number is loaded into TOOSC, the orbital plane at the 
initial time ( t = 0) is used . 

To revoke: 

26015 ELPT OCT 076100000000 NOP 


** Floating point numbers 
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ZMP 


Z 



-XMP 


X 


Figure 3. Instantaneous Lunar Plane Coordinates 


R. Moon Fixed and Rotating Coordinate System 

Vehicle position is provided in two coordinate systems, one 
rotating and one space fixed, based on the earth-moon plane. For 
both systems, the origin is at the center of the earth with the posi- 
tive x-axis directed toward the center of the moon. The x y plane 
is the plane of the moon's orbit about the earth, and the z axis is 
in a northerly direction so as to form a right-hand system. For 
the rotating system (XROT, YROT, ZROT) the x-axis moves with 
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the moon; for the space fixed system (SINJ, YINJ, ZINJ) the x-axis 
is directed toward the position of the moon’s center at injection. 

To use, insert the following modification card: 

25626 BONRCO OCT 076100000000 NOP 

To revoke: 

25626 DONE CO OCT 002000025707 TRA FOSCL 


S . To Print the Coordinates of the Earth Sub-Satellite point and right 
ascension and declination in other than earth or moon reference, the 
following card must be added: 

Modification card: 

14476 NESS OCT 076100000000 NOP 

To revoke: 

14476 NESS OCT 010000014723 TZE SUBSEX 

T. Floating Point Spill Print 

A routine is included which will make the proper correction for 
floating point spills . Normally this routine will print as a warning 

OVERFLOW IN LOC . . . . . , 

To eliminate print, use the following card: 

25006 FPSP OCT 002000025010 TRA *+2 

To restore: 

25006 FPSP OCT 007400403747 TSX E$WOT, 4 

Overflow Diagnostic 

This optional debugging feature prints at each overflow the contents of 
the index registers, AC, MQ, the address part of the three instructions 
preceding the overflow and the numbers in the listed locations. 
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To activate, use the following card: 


25010 OFLDG OCT 076100000000 NOP 

To revoke: 

25010 OFLDG OCT 002000025052 TRA LLXT 


U. Error Return 

If an error is detected the program has three options: 

1 . ERPRT The octal location of the error is written on line and . 

off line. The normal end dump is written on A3 and 
the program goes to the next case. 

2 . ERPRTS The program behaves as above except that it requests 

a memory dump and halts . 

3 . ERPRTC The program behaves as in ERPRT except that it will 

continue the current case through three of these errors 
before going on to next case. 


V . The Flight Path Angle and Azimuth referred to both the geodetic and 

geocentric horizontal planes in earth reference are printed normally. 

To eliminate: 


14722 EAZCAL 
To revoke: 

14722 EAZCAL 

In Moon Reference 
To eliminate: 

14762 MAZCAL 
To revoke: 

14762 MAZCAL 


OCT 076100000000 


OCT 007400436120 


OCT 076100000000 


OCT 007400436120 


NOP 


TSX AZZCAL, 4 


NOP 


TSX AZZCAL, 4 
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W . Additional Stopping Conditions were introduced for moon satellites 
when the radius from the moon is too large, or too small. These 
distances may be changed by the following cards: 


Modification card: 



30226 MRVEC 

DEC 

Minimum distance** 



(.27 ER currently) 

30227 MRVEC+1 

DEC 

Maximum distance** 


(10. ER currently) 


X. Shadow Calculation 


The time boundaries defining the umbra and penumbra are determined 
(including refraction) by a linear interpolation across the two nearest 
integration steps . The radiation pressure switch is controlled only to 
the first integration step after a change occurs. Subsequent updates 
will include an optional refinement to integrate exactly to the boundary. 

The entry, exit, time in and accumulated time for sunlight, umbra and 
penumbra are printed at each change. The following card activates 


shadow testing and printing: 



15334 SHAENT 

OCT 

007400426171 

TSX WSHADE, 4 

To revoke: 




15334 SHAENT 

OCT 

076100000000 

NOP 


Y. Polar Scale Factor is used to change the units of input quantities (for 


polar input only) . 

The modification cards are: 

30145 ICONT DEC . . . . . . (1)** 

30146 POLSC DEC (2)** 


(1) Number of time units in one hour . 

(2) Number of distance units in one earth radius . 

** Floating point numbers 
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To revoke: 


30145 ICONT DEC 3600. , 6378. 165 ** 


Z. Ecliptic Coordinate Output 

This gives additional output (Section X-B-ll) referenced to the ecliptic 
plane. 

To initiate, include: 

26016 ECELPT OCT 002000034740 TRA ECLCOO 

To revoke: 

26016 ECELPT OCT 076100000000 NOP 

A-l . Impact Parameter Calculation (Appendix U) 

This calculates miss distances on the impact plane. 

To use, include: 

12021 TIMPPS OCT 076100000000 NOP 

To revoke: 

12021 TIMPPS OCT 002000012074 TRA ETIMPP 

If several cases are stacked, the same impact plane is used for all of 
them . If a new impact plane is desired, include the card: 

35660 NOMTR DEC 0 

B-l. Geocentric Velocity and Geodetic Position Input 
To use, include: 

12310 LBASC-4 OCT 007400412317 

To revoke: 

12310 LBASC-4 OCT 007400412320 


TSX BASCOR, 4 
TSX POLCAR, 4 


** Floating Point Numbers 
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0-1. Trajectory Search 


The search program is used in connection with the polar load and is 
activated in two ways: 


1 . If the search consists of circumlunar trajectories with 
return to earth 


15204 ITENT 
44001 HELCO 


OCT 002000041307 
DEC x, x, x, x 


TRA STIT 
x = l or 0* 


DEC 1, 0, 0, 0 

DEC 0, 1, 0, 0 

DEC 1, 1, 0 , 0 

DEC 0, 0, 1, 0 

DEC 0, 0, 0, 1 

DEC 0, 0, 1, 1 

etc. 


Outgoing miss parameter only 

Return miss parameter only 

Both miss parameters 

Outgoing osculating elements near perigee 

Return osculating elements near perigee 

Both osculating elements 


2 . If the search is planetary 

15205 VTTENT OCT 002000042334 TRA VIT 

The search input is also loaded through modification cards, a 
sample of which follows: 

43224 MAXIT DEC 5 (maximum number of 

iterations desired)* 

42603 IVAR DEC 0, 0, 10., .0001, 0, 0, 0** 

Input quantities to be varied in the units and sequence of the polar 
load, i.e. , latitude, longitude, altitude, velocity, azimuth, flight 
path angle and time, and the amounts by which to vary the inputs 
to get the partial derivatives are loaded. In the sample given, vary 
height 10 km, speed ,0001 km/sec, and keep all others constant. 

Next 

42614 LQUANT DEC 5, 10* 

the code of the quantities to be achieved. In the sample, pericynthion 
and perigee . The quantity code is as follows: 


* Fixed point numbers 
** Floating point numbers 
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1 . 

Lunar inclination 

(degrees) 

2. 

Lunar ascending node 

(degrees) 

3. 

Lunar argument of pericynthion 

(degrees) 

4. 

Lunar time of pericynthion 

(hours) 

5. 

Lunar pericynthion radius 

(km) 

6. 

Earth inclination 

(or planet) 

7. 

Earth ascending node 

(or planet) 

8. 

Earth argument of perigee 

(or planet) 

9. 

Earth time of perigee 

(or planet) 

10. 

Earth perigee radius 

(or planet) 

11. 

B’T - Miss parameter 

(planet or moon) 

12. 

B’R - Miss parameter 

(planet or moon) 

Next 



42626 QUANT DEC 6000., 5000.** 

are the desired values of the dependent variables . These must 
be listed in the same order as LQUANT. Therefore, these are 
pericynthion and perigee radius (in km), respectively. 

Finally, the tolerances on the above quantities should be given 
and again in the same order as the preceding two, i.e. , 

LQUANT and QUANT 

42640 ICONV DEC 10., 100.** 

also in km. Thus the tolerance on pericynthion radius is 10 km 
and for perigee 100 km. If the solution converges to within the 
specified tolerance, the iterations will stop. 

A maximum of any seven dependent variables can be selected. 
All of the search input modifications must be repeated for each 
repeat or stacked case with appropriate changes as desired. 

To repeat — the sample input for search is as follows: 


15204 

OCT 

002000041307 

ACTIVATE SEARCH 

44001 

DEC 

1, 1, 1, 1 

IIELCO* 

43224 

DEC 

5 

MAXIT* 


* Fixed point numbers 
** Floating point numbers 
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42603 

DEC 

0, 0, 10., .0001, 0, 0, 0 

IVAR** 

42614 

DEC 

5, 10 

LQUANT* 

42626 

DEC 

6000., 5000. 

QUANT** 

42640 

DEC 

10., 100. 

ICONV** 


The normal input for the trajectory search is by modification 
cards . The lunar search option normally iterates on osculating 
elements referred to the equatorial system. To iterate on 
elements referred to the earth-moon plane, use the following 
card: 

41360 EMP6 OCT 076100000000 NOP 

To revoke: 

41360 EMP6 OCT 002000041404 TRA EQCOO 

D-l. Ephemeris Time 

The planetary coordinates are interpolated using ephemeris time 

ET = UT + AT 

An approximate value of AT (35 sec) is used. To change, use the 


following card, giving AT in hours . 

30164 ETMUT DEC (AT inhrs)** 

To revoke: 

30164 ETMUT DEC ,009888888889 (AT is 35 sec) 

E-l. If it is desirable to include additional descriptive material in the title 
printout, the following locations may be used: 

21070 TTPRN BCD Descriptive 

information 

21102 BCD Descriptive 

„ , , , . information 


# _ , . xmwj.mM.wv 

20 symbols (2 cards) maximum. 

The first letter of the first card should be in column 14. 
The first letter of the second card may be in column 13. 


* Fixed point numbers 

** Floating point numbers 



F-l. Beta Integration 


The program has an option for integrating Beta instead of time . Com- 
putation is much faster in this mode but printout doesn't occur at exact 
times . Integration intervals have to be chosen with care . 

Modification card: For Beta Mode 

40537 BETRIG OCT 076100000000 NOP 

To revoke: 

40537 BETRIG OCT 002000000000 TRA 

G-l . Apogee, Perigee and Nodal Crossing Prints are activated by the 
following modification cards: 

30342 CCNT DEC n * 

44214 CANT DEC n * 

44215 CP NT DEC n * 

til 

CCNT prints every n nodal crossing 
CANT controls apogee prints 
CPNT controls perigee prints 


vm-21 



IX. 


INPUT 


This section contains the following information: 

A. General inputs necessary for running a case. 

B . Stacking cases . 

C. Sample input for polar coordinates, no modifications. 

D. Sample input for Cartesian coordinates, no modifications. 

E . Sample input with modifications . 


F. Radar input information. 

A. General Inputs 

Title and Units See Note 

Modification cards, if any (1) 

TRA 2, 4 (2) 

DEC 0, 1, 2, 3 or 4 (3)* 


0 - Cartesian coordinate input 

1 - Equatorial polar coordinate input 

2 - Translunar plane input 

3 - Ecliptic polar coordinate input 

4 - Osculating element input 

DEC 4, 5 or 6 (4)* 

4 - Venus 

5 - Mars 

6 - Jupiter 

are tested for possible reference planets during the flight 
and one must be designated. Moon, Sun or Earth are tested 
automatically and need not be designated. 


* Fixed point numbers. 
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Load 

Number 

Cartesian 

Coordinates 

Polar 

Coordinates 

Translunar 

Plane 

Element 

Load 

See 

Note 

-1 

Omit 

Latitude - 
Degrees and 
decimals 

Lead angle ^ - 
Degrees and 
decimals 

1, 2, 3 or 4* 

(3)** 

0 


Case number 



(5)* 

X 

x -km 

Longitude - 
Degrees and 
decimals 

Argument of 
radius - 
Degrees and 
decimals 

Argument of 
perigee 

(3)** 

2 

y-km 

Altitude - km 

Altitude - km 

RA of ascending 
node 

(3)** 

3 

z -km 

Speed - km/sec 

Speed - km/ sec 

Inclination 

(3)** 

4 

x- km/sec 

Flight path azi- 
muth - Degrees 
and decimals 

Inclination of 
translunar 
plane - Degrees 
and decimals 

Semi-major axis 

(3)** 

5 

y-km/sec 

Flight path angle- 

• Flight path 

Eccentricity 

(3)** 


Degrees and angle - Degrees 

decimals and decimals 


6 z - km/ sec Longitude of Immaterial 

vernal equinox - 
Degrees and 
decimals or 
indicator 


Uses number (3)(15) 

from (-1) 

X - Time of perigee 

2 - Mean Anomaly 

3 - Eccentric Anomaly 

4 - True Anomaly 


* Fixed point numbers 

** Floating point numbers 



Load 


imber 

Title and Units 

See Note 

7 

Launch time - days 

( 8)(6)** 

8 

Launch time - hours 

(8)** 

9 

Launch time - minutes 

(8)** 

10 

Launch time — seconds 

(8)** 

11 

Initial print suppress - hours 


12 

Year of launch 

(17)* 

13 

Initial reference - origin indicator 

1 - Earth 4 - Venus 

2 - Moon 5 - Mars 

3 - Sun 6 - Jupiter 

<3)<9)* 

14 

Thrust indicator 

0 - No thrust 

1 - Thrust 

(10)* 

15 

Maximum time of flight - hours 

(11)** 


If thrust is used, the number 1 has to be loaded into M Load+14" and 
additional input data have to be supplied. These cards are inserted following 
the regular input, but preceding the radar input, if any . The thrust load has 


the following format: 

TO Integration interval to be used during thrust period (16)** 

T1 Number of thrust periods, 24 maximum (18)* 

T2 Initial mass of vehicle - pounds (19)** 

T3 Mass flow-pound mass/ pound thrust/hour ** 

T4 Initial time of first thrust period *.* 


♦Fixed point numbers. 
♦♦Floating point numbers. 
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For each thrust period: 


1T1 

T 

X 

x-component of thrust - pounds 

** 

1T2 

T 

y 

y-component of thrust - poinds 

** 

1T3 

T 

z 

z-component of thrust - pounds 

** 

1T4 

0 , 0 

Space reserve 


1T5 


Termination time of this thrust period 

** 


If simulated radar indormation is desired, the following input must be 


supplied: (See VTH-C) (12) 

RO Number of stations (13)* 

R1 BCD 4 - Station name (5)(14) 

R2 Longitude - degrees, minutes, seconds (12)(14)** 

R3 Latitude - degrees, minutes, seconds (12)(14)** 

R4 Altitude - feet (12)(14)** 


NOTES : 

1. See Section Vm for modifications available. 

2 . Terminates reading of modification cards . This card must always be 

included. 

3. See Sections IV and VUI and translunar plane section. 

4. The choice of planet is significant only if trajectories are to be 

computed which will come close to either Mars or Venus or Jupiter. 

In this case, the planet designated may serve as a reference body 
if its sphere of influence is entered. For all other trajectories, the 
designation of a planet only affects the output and is usually unimportant. 
A choice, however, must be made. 

* Fixed point number 
** Floating point numbers 
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NOTES : (cont.) 


5. These inputs are used only for identification and not for computing 

purposes . 

6. Number of days since December 31, 0. 0 hrs UT. (See footnote 1 of 

Section IV.) 

7. This suppresses printing until after specified time has been exceeded. 

The input time is arbitrary and may be chosen as desired . 

8 . The number of days must be a floating point integer, the minutes and 

seconds may be included as a decimal fraction of hours . In this case, 

0 must be entered on cards 9 and 10. 

9. Indicates the body whose center is the origin of the initial reference 

system. 

10. For a more detailed description of the use of the thrust programs, see 

Section IX- A, B (Input) . 

11. See Section V. 

12. See Section VH. 

13. Enter number of tracking stations . 

14. Enter the following information for each station: (See Section VII) 

BCD 4 Station name, 24 symbols maximum 
DEC Longitude 

DEC Latitude 

DEC Altitude 

Longitude and latitude may be given in degrees and decimals if more 
convenient. 0 then has to be used in place of minutes and seconds. 

The resulting printout of the polar station coordinates will then show 
the integral number of degrees only . The station position in core, 
however, will reflect the fractional part. 

15. For polar coordinate input only. If zero is loaded in this position, the 

load in position 1 'longitude" is interpreted as right ascension. If 
a fixed point 1 is loaded, the program computes the longitude of the 
vernal equinox at launch time from the data supplied on cards 7, 8, 9 
and 10. If for some reason this computation is not desired, a floating 
point number is loaded on card 6 and then this number will be used as 
the longitude of the vernal equinox. 
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NOTES : (cont. ) 


16. If fractional parts of hours are used for integration intervals, it is 

recommended that multiples of negative powers of 2 are used to 
eliminate round-off error in the time. 

17. This quantity ensures that the proper file on the planetary tape is used. 

It will stop the program if the wrong tape is mounted or the proper 
file cannot be found. 

18. Periods of no thrust intervening between thrust periods must be in- 

cluded in this count. 

19. No staging is considered, but may be included if required. 

20. For a period of no thrust intervening between two thrust periods, 

cards 1T1 through 1T4 have to be included with 0 in place of 

T , T and T . In each case the final time t, is used as initial 
x v z f 

time for the next thrust period, until the indicated number of thrust 

periods has been reached. 


B. Stacking of Cases 

If many cases are to be run differing only by a few parameters, it is 
not necessary to repeat that part of the input which remains unchanged. The 
changes only are to be loaded but an octal location must be punched in columns 
2 through 6. For instance, if only the latitude is to be changed, the following 
card is loaded: 

31263 DEC ...... latitude 

The table on the following page gives a complete list of the load and their octal 
locations . 

If the launch time is changed in a polar load stack, the computation of 
the longitude of the vernal equinox must be re-initiated by loading 1 into LOAD+6 . 

These cards are followed by three TRA 2, 4 to terminate loading. If 
radar is used and the station coordinates are unchanged, an "Avoid Radar Load" 
card (described in Section Vm-C) has to be included. If thrust is used, the 
thrust load has to be repeated or the following modification card has to be in- 
cluded to cancel the thrust loading: 
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TABLE 1 


LOAD LOCATIONS 


Load 

Number 

Octal 

Location 

Title 

- 1 

31263 

latitude 

0 

4 

case number 

1 

5 

x or longitude 

2 

6 

y or height 

3 

7 

z or speed 

4 

31270 

x or azimuth 

5 

1 

y or flight path angle 

6 

2 

z or longitude of vernal equinox 

7 

3 

days 

8 

4 

hours 

9 

5 

minutes 

10 

6 

seconds 

11 

7 

print suppress 

12 

31300 

year 

13 

1 

origin indicator 

14 

2 

thrust indicator 

15 

3 

maximum time 
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12630 THRST+2 


OCT 076100000000 


NOP 


To revoke: 


12630 THRST+2 


OCT 007400420407 


TSX THLD, 4 


Sample Input for Polar Coordinates - with repeat case 


30153 DEC 1 


TESSTR - modification to include 
tess.er.al harmonics 


21070 BCD 


FRED SHAFFER 


TRA2.4 


Go to regular load - end of modification 
cards 


DEC 1 


Polar coordinates 


DEC 5 


Mars is planet choice 


DEC 30. 60860 Latitude in degrees and decimals 


DEC 1 


Case number 


DEC - 71.78139 Longitude in degrees and decimals 


DEC 160. 72594 Altitude in kilometers 


DEC 7.8439280 Velocity in km/sec 


DEC 77.89589 


Flight path azimuth in degrees and 
decimals 


DEC - .0317 


Flight path angle in degrees and 
decimals 


DEC 1 


Program computes longitude of vernal 
equinox at launch time 


DEC 268. 


Launch time in days from Dec. 31, 
0, Oh UT 
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DEC 

14. 

Launch time in hours from midnight UT 

8 


DEC 

0 . 

Launch time in minutes 

9 


DEC 

0 . 

Launch time in seconds 

10 


DEC 

0 . 

Print suppress in hours 

11 


DEC 

1962 

Year of launch 

12 


DEC 

1 

Earth is initial reference origin 

13 


DEC 

0 

No thrust 

14 


DEC 

24. 

Maximum time of flight in hours 

15 

30153 

DEC 

0 

No tesseral harmonics in repeat case 


21102 

BCD 


NO TESSERAL HARMONICS 



TRA 

2,4 




TRA 

2,4 




TRA 

2,4 




TRA 

37777 

End of run, non-monitor 



D. Sample Input for Cartesian Coordinates - no modifications 


TRA 2,4 



DEC 0 

Cartesian coordinates 


DEC 4 

Venus is planet choice 


DEC 1 

Case number 

0 

DEC 4953.3445 

x-coordinate in km 

1 

DEC 3716.7623 

y-coordinate in km 

2 

DEC 2224.0628 

z -coordinate in km 

3 
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DEC 

7.1423107 

x-component in km/ sec 

• 4 

DEC 

7.6362243 

y- component in km/sec 

5 

DEC 

3.1457277 

z -component in km/sec 

6 

DEC 

33. 

Launch time in days 

7 

DEC 

2. 

Launch time in hours 

8 

DEC 

59. 

Launch time in minutes 

9 

DEC 

42.72 

Launch time in seconds 

10 

DEC 

0 . 

Print suppress in hours 

11 

DEC 

1966 

Year of launch 

12 

DEC 

1 

Earth is initial reference origin 

13 

DEC 

0 

No thrust 

14 

DEC 

168. 

Maximum time of flight in hours 

15 

TRA 

r> 

CO 

End of data deck, non-monitor 



E. 


Sample Input for Cartesian Coordinates - with modifications 


12654 

OCT 

007400423630 

Load radar data 

14771 

OCT 

007400424002 

Compute radar output 

25006 

OCT 

002000025010 

Eliminate overflow print 

21070 

BCD 


IMP ORBIT WITH RADAR 

30177 

DEC 

4.5 

Change near, far earth boundary 


TRA 

2,4 



DEC 

0 

Cartesian coordinates 


DEC 

5 

Mars is planet choice 
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DEC 

102 

Case number 

0 

DEC 

6336.4095 

x-coordinate in km 

1 

DEC 

812.60758 

y-coordinate in km 

2 

DEC 

1600. 4426 

z -coordinate in km 

3 

DEC 

.12118217 

x-component in km/sec 

4 

DEC 

9.4967498 

y-component in km/sec 

5 

DEC 

-5.3016564 

z -component in km/ sec 

6 

DEC 

152. 

Launch time in days 

7 

DEC 

11. 

Launch time in hours 

8 

DEC 

3. 

Launch time in minutes 

9 

DEC 

58. 

Launch time in seconds 

10 

DEC 

0 

Print suppress in hours 

11 

DEC 

1963 

Year of launch 

12 

DEC 

1 

Earth is initial reference origin 

13 

DEC 

0 

No thrust 

14 

DEC 

160. 

Maximum time of flight in hours 

15 

Radar Input Information 



DEC 

1 

Number of stations 

R0 

BCD 


4 MID ATLANTIC (station name) 

R1 

DEC 

-30., 0., 0. 

Longitude of station - degrees, 
minutes, seconds 

R2 

DEC 

5., 0., 0. 

Latitude of station - degrees, 
minutes, seconds 

. R3 

DEC 

0.0 

Altitude of station - feet 

R4 

TRA 37777 

End of data deck, non-monitor 
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X. 


OUTPUT 


A. Program Outputs 

The following information is printed as the output of the program. 

1 . Modification cards 

2. Title 

3. Case number and any identifying titles 

4. Launch time - days, hours, minutes, seconds 

5. Input data - printed out exactly as they entered the program 

6 . Octal locations of LOAD-1 through LOAD+7 

7 . List of parameters used in run 

8 . At each rectification the following data are printed: 

RECTIFICATION PRINT (a) REFERENCE 

(b) PERT OVER UNPERT = (c) 

AT TIME = (d) DELTA T = (e) 

(a) Reference body 

(b) and (c) indicate the reason for rectification: 

If (c) - 0 then reference body has been changed as indicated 
by (b) and the following nomenclature is used to indicate the 
new reference body: 

MR - Moon reference VR - Venus reference 

SR - Sun reference MA - Mars reference 

ER - Earth reference JR - Jupiter reference 

If (c) ± 0 then the position, velocity, acceleration perturba- 
tions or the incremental eccentric anomaly have exceeded the 
permissible limits and (b) indicates which has been exceeded 
(see Section VI). These indications are given as: 

PO - Position 

VL - Velocity 

AC - Acceleration 

TH - Incremental eccentric anomaly 

(d) Rectification time in hours from time of launch 

(e) Integration interval to be used 
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9 . At every print time the following data are printed: 

TIME IN DAYS, HRS, MINS, SECS, 

T = 

XR YR ZR RR ______ 

XRDT YRDT ZRDT RRDT 

RIGHT ASCENSION (DEG) = DECL = 

EARTH SUBSAT POINT LONG ~ ______ 

LAT = 

HT = 

GHA = 

(a) Days, hours, minutes, seconds from time of launch 

(b) Time in hours from time of launch 

(c) Position coordinates and magnitude of radius vector 

with respect to the reference body - km 

(d) Velocity components and magnitude of velocity vector 

with respect to the reference body - km/ sec 

(e) Right ascension and declination in earth reference 

system - deg 

(f) Longitude or sub -satellite point - deg 

(g) Latitude (geodetic) - deg 

(h) Geodetic height above the earth’s surface - km 

(i) Greenwich hour angle - deg 

MOON SUBSAT POINT LONG - 

LAT = 

AZIMUTH = 

ELEVATION - ____ • 



OSCULATING ELEMENTS AT TIME T = 


TRUE ANOMALY = 

(n) 

SEM MAJ AXIS - 

(o) 

ECCENT 

(P) 

PERICENT 

(q) 

APOCENT 

(r) 

INCLINATION 

(s) 


(j) Moon longitude - Angle between the projection of the 

vector from the moon to the vehicle onto the moon's 
orbital plane and the moon-earth vector (moon 
reference only) - deg 

(k) Moon latitude - Angle between the radius vector con- 

necting the moon and the vehicle and its projection 
onto the orbital plane of the moon about the earth 
(moon reference only) - deg 

(-t) Selenocentric flight path azimuth - deg 

(m) Selenocentric flight path angle - deg 

(n) True anomaly - deg 

(o) Semi-major axis of trajectory - ER 

positive = ellipse 
negative = hyperbola 

(p) Eccentricity of trajectory** 

(q) Closest distance to the reference body (not necessarily 

the earth) - km** 

(r) Farthest distance from the reference body (not neces- 

sarily the earth) - km** (Meaningful only for elliptic 
orbits . ) 

(s) Inclination of the orbital plane defined as the angle be- 

tween the positive polar axis and the angular momen- 
tum vector - deg** 


** These are the osculating values and hence only constitute an estimate of 
the quantities described. 


X-3 



AEG PERIC 


(t)' 

PERIOD = (u) 

MEAN MOT = ___ (v) 

R A ASC NODE = (w) 

M ANOMALY = (x) 

E ANOMALY = (y) 

T PERIC “ (z) 

UNIT PERICENTER POSITION VECTOR = (aa) 

UNIT ANGULAR MOMENTUM VECTOR = (bb) 


(t) Argument of pericenter - Angle measured from the ascending 

node to the pericenter vector - deg** 

Set to zero for circular orbits and poorly determined for 
near-circular orbits 

(u) Period - hrs** 

(v) Mean motion - rad/hr** 

(w) Right ascension of the ascending node measured from the vernal 

equinox eastward along the equator - deg** 

(x) Mean anomaly - rad** 

(y) Eccentric anomaly - rad** 

(z) Time of nearest pericenter - hrs** 

(aa) Components of the unit vector directed from reference toward 
pericenter** 

(bb) Components of the unit angular momentum vector 


** Osculating values 
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B. 


Optional Outputs 


XVE = 

YVE = 

ZVE = 

RVE = 

(a) 

XVM = 

YVM = 

ZVM = 

RVM = 

(b) 

XME - 

YME = 

ZME - 

RME = 

(c) 

XVS = 

YVS = 

ZVS = 

RVS = 

(d) 

XWN = 

YWN = 

ZWN = 

RWN = 

(e) 

XVMR = 

YVMR = 

ZVMR = 

RVMR = 

(f) 

XV JP = 

YVJP = 

ZVJP = 

RVJP = 

(g) 

XI = 

ETA = 

ZETA = 

PERT = 

(h) 

XEDT = 

ETADT = 

_ ZETADT= 

VPERT = 

(i) 

D2X3 = 

D2ETA = 

_ D2ZETA = 

APERT = 

(j) 


The above optional output appears between XRDT and right ascension 
in the standard output. For instructions on how to obtain, see Section 

vm-o. 

(a) Coordinates of vehicle with respect to the earth - km 

(b) Coordinates of vehicle with respect to the moon - km 

(c) Coordinates of moon with respect to the earth - km 

(d) Coordinates of vehicle with respect to the sun - km 

(e) Coordinates of vehicle with respect to Venus - km 

(f) Coordinates of vehicle with respect to Mars - km 

(g) Coordinates of vehicle with respect to Jupiter - km 

(h) Perturbation vector and magnitude of the perturbations with 

respect to the reference body - km 

(i) Perturbation velocity vector and magnitudes - km/sec 

(j ) Perturbation acceleration vector and magnitudes - km/sec 2 


X-5 



2. Moon Rotating and Fixed Coordinate Systems 


XDSfJ - 

YINJ = 

ZINJ = 

RINJ = 

(a) 

XROT = 

YROT = 

_ ZROT = 

RROT - 

_ (b) 

XEM = 

YEM = 

ZEM - 

REM = 

__ (c) 


(a) Coordinates of the vehicle with respect to a fixed coordinate 

system defined at injection (see Appendix Q) - km 

(b) Coordinates of the vehicle with respect to a rotating coordinate 

system (see Appendix Q) - km 

(c) Coordinates of the earth with respect to the moon - km 


3. Moon Osculating Elements 

REFERENCE PLANE IS MOON ORBITAL PLANE AT (a) 

TIME = 

INCLINATION =_ (b) 

ASC NODE = (c) 

ARG PERICENT = (d) 

CHNG ASC NODE = _____ (e) 

INCL MOON PLANE = ______ (f) 

RA ASC NODE OF MOON = (g) 

ARG OF LAT OF MOON =* (h) 


(a) Time in hours from launch. The reference plane is the oscu- 

lating orbital plane of the moon . The x-axis is along the 
moon-earth vector . The reference system is variable or 
fixed, depending on the option selected (see Section VIII -Q) . 

(b) Inclination of vehicle plane with respect to the orbital plane 

of the moon - deg 
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(c) Ascending node - Angle between the nodal line and x-axis - deg 

(d) Argument of peric enter - deg 

(e) Change in ascending node from time zero - deg 

(f) Inclination of the moon orbital plane with respect to the earth's 

equatorial plane - deg 

(g) Right ascension of the ascending node of the moon - deg 

(h) Argument of latitude of the moon - angle between the moon's 

position vector and the ascending node of the moon - deg 

4. Angle from Ascending Node to Vehicle 

ANGLE FROM ASCEND NODE TO SAT - (a) 

(a) Angle from the ascending node to the satellite - deg 

5. Error Print 

ERROR IN LOC _________ (a) 

(a) Octal location subsequent to where an error occurs . 

6. Floating Point Spill Print 

OVERFLOW IN LOC ____________ (a) 

(a) Octal location of the instruction which caused a floating point 
spill to occur. 

7. End Print 

Four blocks of storage that are useful in checking and debugging. 
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8 . Shadow Print 


PASSAGE FROM 


' SHADOW 


'i 

PENUMBRA 

1 PENUMBRA 

> TO < 

SHADOW 

] SUN 

PENUMBRA 

PENUMBRA 
k > 


SUN 


AT (a) TIME IN < 


SHADOW 'I 
PENUMBRA l 
SUN 

PENUMBRA 


(b) ACCUMULATED TIME (c) , (d) 


(a) Time at which vehicle crosses denoted shadow boundary - hrs 

(b) Total time the vehicle spends in denoted shadow region during 

current traverse - hrs 

(c) Total accumulated time spent in denoted shadow region since 

launch - hrs 

(d) 0 indicates earth’s shadow 
1 indicates moon's shadow 


9. Radar Output 

STATION 

AZIMUTH 

ELEVATION 

TOPOC. RA 
TOPOC. DECL. 

SLT RNG 

RANGE 


(a) 

(b) 

(b) 

(c) 

(c) 

(d) 

(e) 


(a) Station name (identification) for each station 

(b) Azimuth and elevation with respect to each station - deg 

(c) Topocentric right ascension and declination with respect to each 

station - deg 

(d) Slant range to each station - km 

(e) Rate of change of slant range for each station - km/ sec 

If the elevation is negative (vehicle is below horizon), this 
print is suppressed for the station of reference. 
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10. Flight Path Azimuth and Angle Output 


GEOCENTRIC AZIMUTH = _____ (a) 

ELEVATION = (b) 

GEOD AZIM = (c) 

ELEV = (d) 


(a) Geocentric flight path azimuth - deg 

(b) Geocentric flight path angle - deg 

(c) Geodetic flight path azimuth - deg 

(d) Geodetic flight path angle - deg 

11, Ecliptic Coordinate Output 


EXES = 

EYES = 

EZES 

ERES = (a) 

ELONG = Jbl 

EXTS * 

EYTS * 

EZTS =_ 

ERTS = (c) 

. TLONG = Jdi 

EXVS = 

EYVS - 

EZVS=_ 

ERVS = (e) 

. VLONG = Jf)_ 

XECL =_ 

YECL = 

ZECL = 

RECL= (g) 


DXECL = 

DYECL == 

DZECL = 

DRECL = (h) 



(a) Ecliptic coordinates of earth with respect to the sun - km 

(b) Ecliptic longitude of earth - deg 

(c) Ecliptic coordinates of planet with respect to the sun - km 

(d) Ecliptic longitude of planet - deg 

(e) Ecliptic coordinates of vehicle with respect to the sun - km 

(f) Ecliptic longitude of vehicle - deg 

(g) Ecliptic coordinates of vehicle with respect to reference body - km 

(h) Ecliptic velocity of vehicle with respect to reference body - km/ sec 

ANGLE BETWEEN RVE-RVS = til (j) 

ELEMENTS REFERRED TO ECLIPTIC PLANE AT TIME T = 

ECC = (k) 
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INCL = . (1) 

ASCN - (m) 

ARPG = (n) 

LONGRREF = (o) 

LATRREF = ___ (p) 

AZIMUTH = (q) 

ELEV = (r) 

( i ) Angle between vehicle -earth and vehicle -sun vectors - deg 

( j ) Angle between vehicle -planet and vehicle -sun vectors - deg 

(k) Eccentricity of trajectory 

(l) Inclination of vehicle plane with respect to the ecliptic plane - deg 

(m) Ascending node - angle between nodal line and x-axis - deg 

(n) Argument of perigee - deg 

(o) Ecliptic longitude of vehicle with respect to the reference body - deg 

(p) Ecliptic latitude of vehicle with respect to the reference body - deg 

(q) Flight path azimuth with respect to the reference body on celestial 

sphere - deg 

(r) Flight path angle with respect to the reference body on celestial 

sphere - deg 


12. Reentry Output 

REENTRY PRINT TIME INERTIAL SPEED (km/sec) 

Right ascension, declination, earth subsatellite points, and flight 
path azimuth and angle as given above, 

13. Trajectory Search Output 

The output consists of the normal ITEM output for a nominal 
trajectory and the same trajectory output for each variation re- 
quested. The output format used only for the trajectory search 
follows: 
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VARIATION IN INITIAL 

CONDITIONS (a) (b) (c) (d) 

QUANTITY CODE 

DESIRED VALUES OF ABOVE QUANTITIES 
REQUIRED ACCURACY 

(a) Change in latitude - deg 

(b) Change in longitude - deg 

(c) Change in altitude - km 

(d) Change in velocity - km/sec 

(e) Change in azimuth - deg 

( f ) Change in flight path angle - deg 

(g) Change in initial time - hrs 

(h) Code indicating quantities (up to 7) to be searched for 

( i ) Desired values of above quantities 

( j ) Tolerances allowed on above values 

MATRIX OF PARTIAL DERIVATIVES (k) 

RESIDUALS AND CHANGES IN INITIAL 

CONDITIONS q> (m) (n) 

(k) Matrix having the dependent variables arranged by row; the 

independent by column 

(•i) Residuals (desired-nominal) of quantities designated by the 
quantity code 

(m) Changes required in initial conditions 

(n) Normalized changes in initial quantities in order of the variations 
14. Impact Parameter Output 

IMPACT PLANE PARAMETERS DIRECTION COSINES OF 


- is ) m ssl 

ilD_ 

jU1_ 

(j) 


T VECTOR 
R VECTOR 


M 

M 



ASYMPTOTE 


(e) (f) 


MISS PARAMETERS CROSS T 
B - T , B-R 

(a) Unit vector parallel to the ecliptic or moon orbital plane and per- 

pendicular to the incoming asymptote 

(b) Unit vector perpendicular to T-vector and the asymptote in a 

right-hand sense 

(c) Unit vector in the direction of the incoming asymptote 

(d) Time at which the impact plane is crossed - hrs 

(e) The dot product of R^.^, at crossing time and T-vector - km 

( f ) The dot product of Ry^ at crossing time and R -vector - km 

(a), (b) and (c) are printed after the set up of the impact plane 
coordinate system. This occurs when the nominal trajectory reaches 
the required approach distance. This will not be repeated for sub- 
sequent trajectories unless NOMTR is reset to zero. 

(d), (e) and (f) are printed after the spacecraft crosses the 
impact plane . 

15. Overflow Diagnostic 

ADDRESS PART OF 3 PRECEDING INSTRUCTIONS 

AND IR STATUS (a) 

_ib}_ 

NUMBERS IN LISTED LOCATIONS, 

AC, MQ (c) (d) (e) 

(a) Address part of the three instructions preceding the location of 

the overflow 

(b) The contents of index registers 1, 2, 4 at the time of the overflow 

(c) The three numbers in the locations referred to in (a) 

(d) The contents of the accumulator 

(e) The contents of the multiplier quotient register 
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16.. Earth-Fixed Velocity 


EARTH-FIXED VELOCITY X = (a) Y = (b) Z = (c) 

V - (d) AZ = (e) EL = (f) 

(a) x-component of the earth-fixed velocity - km/sec 

(b) y-component of the earth-fixed velocity - km/sec 

(c) z-component of the earth-fixed velocity - km/sec 

(d) Magnitude of earth-fixed velocity - km/sec 

(e) Geocentric flight path azimuth - deg 

(f ) Geocentric flight path angle - deg 
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XI. . INTERNAL PROCEDURES 
' A. Units 

The units which are used internally are earth radii and earth 
radii per hour in earth and moon references, and astronomical units 
and astronomical units per hour in sun and planet reference systems. 


B. Ephemeris Tape 

1. The relative positions of the solar system bodies are obtained 
from cards furnished by the U. S. Naval Observatory. A separate pro- 
gram prepares a binary tape referred to the mean equinox of date with 
15 days per record in a form compatible with the main program. The 
main program searches the tape and reads in the proper file and record, 
keeping 30 days of tables in core storage at a time. 

The first record on each file consists of the year in fixed decimal. 
Each of the following records contain the following information: 

Word 1 : Initial time of record in hours from base time. 

Then 12 consecutive 15 word blocks (one word per day) containing 
the equatorial coordinates of: 


XSE 

YSE 

ZSE 

Sun with respect to earth 

XJS 

YJS 

ZJS 

Jupiter with respect to the sun 

XAS 

YAS 

ZAS 

Mars with respect to the sun 

XVS 

YVS 

ZVS 

Venus with respect to the sun 

Then three 30 

word blocks 

containing: 

XME 

YME 

ZME 

Moon with respect to earth 


h h 

The moon coordinates are stored in half-day intervals (0.0 ,12 .0 UT). 
The unit of distance is the radius of the earth. All other tables are in 
daily intervals (0* 1 . 0 UT) and the unit of distance is AU. 

At present, an ephemeris tape is available for 1961-1970, written 
in nine, two-year files. The files overlap one year. 
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2. The tape generated from JPL data is written in two-year over- 
lapping files covering the period 1968-1982. The first record of each 
file has the starting year twice, then the number of records and, 
finally, the number of files. Each succeeding record contains the 
following: 

Time in hours with respect to the beginning of the file. 

Mercury coordinates with respect to the sun. 

9 x-values 4 days apart 

9 y-values 4 days apart 

9 z-values 4 days apart 

Venus coordinates with respect to the sun. 

5 x-values 8 days apart 

5 y-values 8 days apart 

5 z-values 8 days apart 

These are followed by the coordinates of the sun with 
respect to the earth, Mars with respect to the sun, 

Jupiter with respect to the sun, Saturn with respect to 
the sun, Uranus with respect to the sun, Neptune with 
respect to the sun, Pluto with respect to the sun, and 
earth-moon barycenter with respect to the sun, all 
eight days apart. 

The positions of the moon with respect to the earth 
follow: 

33 x-values half-day intervals 

33 y-values half-day intervals 

33 z-values half-day intervals 

These position values are then followed by their respec- 
tive velocities in the same order. 


C. Ephemeris in Core 

1. The astronomical tables are stored in core using 24-hour inter- 
vals for the sun and planets, and 12-hour intervals for the moon. 

There are always 30 days of tables available, arranged in such a way 
that the value of time for which the interpolation takes place is not near 
either end of the table. 

In location TABLE is stored the time of the first entry from 
the initial time. Stored in the following locations are: 

TABLE + 1 to TABLE + 30 x coordinates of the sun 
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TABLE + 31 to TABLE + 60 
TABLE + 61 to TABLE + 90 
TABLE + 91 to TABLE + 180 
TABLE + 181 to TABLE + 270 
TABLE + 271 to TABLE + 360 
TABLE + 361 to TABLE + 420 
TABLE + 421 to TABLE + 480 
TABLE + 481 to TABLE + 540 


y coordinates of the sun 
z coordinates of the sun 
x, y , z coordinates of Jupiter 
x,y, z coordinates of Mars 
x, y , z coordinates of Venus 
x coordinates of the moon 
y coordinates of the moon 
z coordinates of the moon 


2, Using the JPL ephemeris tape, the astronomical tables stored 
in core are: 

TABLE time 

TABLE + 1 to TABLE + 9 XSE 

TABLE + 10 to TABLE + 18 YSE 

TABLE + 19 to TABLE + 27 ZSE 

TABLE + 28 to TABLE + 36 XJPS 

TABLE + 37 to TABLE + 45 YJPS 

TABLE + 46 to TABLE + 54 ZJPS 

TABLE + 55 to TABLE + 63 XMRS 

TABLE + 64 to TABLE + 72 YMRS 

TABLE + 73 to TABLE + 81 ZMRS 

TABLE + 82 to TABLE + 90 XVNS 

TABLE + 91 to TABLE + 99 YVNS 

TABLE + 100 to TABLE + 108 ZVNS 

TABLE + 109 to TABLE + 173 XME 

TABLE + 174 to TABLE + 238 YME 

TABLE + 239 to TABLE + 303 ZME 
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The velocities are stored in the same manner as follows: 


TABLE + 304 to TABLE + 312 XDSE 

TABLE + 313 to TABLE + 321 YDSE 

TABLE + 322 to TABLE + 330 ZDSE 

TABLE + 331 to TABLE + 339 XDJPS 

TABLE + 340 to TABLE + 348 YDJPS 

TABLE + 349 to TABLE + 357 ZDJPS 

TABLE + 358 to TABLE + 366 XDMRS 

TABLE + 367 to TABLE + 375 YD MRS 

TABLE + 376 to TABLE + 384 ZDMRS 

TABLE + 385 to TABLE + 393 XDVNS 

TABLE + 394 to TABLE + 402 YDVNS 

TABLE + 403 to TABLE + 411 ZDVNS 

TABLE + 412 to TABLE + 476 XDME 

TABLE + 477 to TABLE + 541 YD ME 

TABLE + 542 to TABLE + 606 ZD ME 


Perturbation Program 


The perturbation program solves three differential equations 
for XI, ETA, ZETA. The differential equation for XI, with the 
various terms replaced by the storages containing them, is repre- 
sentative of all three equations and is given here. 

D2XI = - GME (VCOR+ 0/VCOR + 3 - COMP + 0/COMP + 3) 

- GMVN ( VCOR + 18/VC OR + 21 - COMP + 18/COMP + 21 ) 

- GMS (VCOR + 12/VCOR + 15 - COMP + 12/COMP + 15 ) 

- GMMR (VCOR + 24/VC OR + 27 - COMP + 24/COMP + 27 ) 

- GMJP (VCOR + 30/VCOR + 33 - COMP + 30/COMP + 33 ) 

- GMM (VCOR + 6/VCOR+ 9 - COMP + 6/COMP + 9) 

+ OTHER PERTURBATIONS 
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where, e. g. , in the first term GME = K times the mass of the earth, 
and VC OR + 3 is the length cubed of the vector (VC OR + 0, VC OR + 1, 
VCOR + 2 ). Similarly, in the other terms the denominator is the length 
cubed of the corresponding numerator. In case the two terms within 
each parenthesis are nearly equal, they are computed by the special 
method described in Appendix E to avoid loss of accuracy. The contents 
of the COMP storage at any time, t, depends upon the reference origin 
at that time. For details see Tables 3 and 4. 

Here XME refers to the x-component of the vector from earth 
to moon, with corresponding definitions for the other quantities. An 
additional subscript of 0 to any quantity denotes the vector derived 
from the two-body problem. 
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TABLE 3 - CONTENTS OF COMP TO COMP + 33 
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XII. OPERATING PROCEDURES 

A. Introduction 

The following procedures for running this program are assembled 
to serve as a guide for programmers and IBM 7094 machine operators 
when using the Interplanetary Trajectory Encke Method program. It is 
strongly advised that, before using the program, all input should be listed 
and examined carefully , then fastened to output when received in order to 
provide a permanent record of input for each case run. 


B. Normal Operation (Non-Monitor) 

1. Input ; The data are supplied on cards. They should be 
loaded on tape and mounted on A2. The last card of the 
input data should be a TRA 37777 (for non-monitor 
operation) beginning in column 8. The normal program 
stop in this case will show location 37777 in the instruc- 
tion counter and HTR 37777 in storage. 

2. Output ; The output of the program is on tape A3. 

3. Ephemeris Tapes ; A variety of ephemeris tapes (tables 
of planetary positions as functions of time and year) may 
be used. Normally a 10-year tape will be used. The 
particular tape desired will be noted from the input in- 
structions by tape number. These tapes may be either 
high or low density and should be so marked on the tape 
reel along with the other identifying information. The 
ephemeris tape should be mounted on A5 with the appro- 
priate density. These tapes should all be FILE PROTECTED 


C. Monitor Operation 

For monitor operation, the input card TRA 37777 is omitted and 
the program followed by the input is loaded off-line on A2. The end 
of file mark on A2 after all cards are read will then return control to 
the monitor. The loader (first card) must be replaced by an A2 loader. 

The instructions on the call card are as follows: 
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00000 

0 

00025 

0 

00003 

IOCD 

00001 

0 

06000 

0 

00001 

TCOA 

00002 

0 

76000 

0 

00004 

ENK 

00003 

0 

13100 

0 

00000 

XCA 

00004 

0 

73400 

1 

00000 

PAX 

00005 

0 

77200 

0 

02205 

REWB 

00006 

-2 

00001 

1 

00014 

TNX 

00007 

0 

76200 

0 

02225 

RTBB 

00010 

-0 

54000 

0 

00022 

RCHB 

00011 

0 

06100 

0 

00011 

TCOB 

00012 

2 

00001 

1 

00007 

TIX 

00013 

-0 

03000 

0 

00014 

TEFB 

00014 

-0 

02200 

0 

00015 

TRCB 

00015 

0 

76200 

0 

02225 

RTBB 

00016 

-0 

54000 

0 

00021 

RCHB 

00017 

-0 

54400 

0 

00000 

LCHB 

00020 

0 

02000 

0 

00001 

TRA 

00021 

-1 

00003 

0 

00000 

IOCT 

00022 

2 

00000 

0 

00000 


00023 

1 

00000 

0 

00022 



Operating Procedure 

a. Place input tape on A2 

b. Place output tape on A3 

c. Place ephemeris tape on A5 

d. Place program tape on B5 

e. All sense switches up 

f. Key number 35 down (for first file on B5) 

g. Clear machine 

h. Ready B5 call card in card reader 

i. Press load cards key 

Program will start and continue to normal stop; IC 37777 and 
HTR 37777 in storage. 


5 

14, 1, 1 
5 
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D. Alternate Method of Operation 


1. Occasionally it is useful to place the data immediately behind the 
binary deck and load the program directly into the machine from 
tape A2. The A2 loader must replace the B5 loader as the 
first card in the deck and the A2 call card is used to initiate the 
program. 

2. To run with deck loaded on A2 

a. Place A2 loader on binary deck 

b. Place data behind binary deck 

c. Load cards onto tape; mount on A2 

d. Place planetary ephemeris tape on A5 

e. Key 35 down 

f. Place A2 call card in card reader 

g. Press load cards key 

h. Output is on A3 


E. Error Stops 

Numerous error stops are provided when the machine detects in- 
consistencies. In this case, the machine prints ERROR AT LOCATION 
........... , then prints a number of working storages to aid in 

finding the cause of the error and proceeds to the next case. Most 
frequently , such conditions are caused by faulty input. As a first step, 
therefore, the input should be examined carefully. In this procedure 
it is extremely helpful to list the input before executing the program 
and clipping it to the oulput. 


F. Machine Requirements 

The program uses a 32K, 2 -channel IBM 7090 or 7094 computer. 
It uses 3 tapes and 1 printer on A channel and 1 tape on B chan- 
nel. A second tape is used on B channel if the save-restore feature 
is used. For monitor operation, or the alternate mode described in 
Section XII-D, the B channel is unnecessary. 


G. Additional Information 


For additional information to aid in operating this program, please 
contact Mr. F. H. Whitlock, Goddard Space Flight Center, Theoretical 
Division, Greenbelt, Maryland 20771. 
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XHI. LIST OF IMPORTANT LOCATIONS IN THE PROGRAM 


This section presents a list of locations of constants and other quantities 
which will enable the program user to modify the computations for special 
purposes. 

The majority of the solar system constants, e.g. , earth gravitational 
harmonics, planetary masses, etc. , are found in the program listing beginning 
at 25415. The first group of coefficients are used in the series expansion for 
the sub-satellite computation. Notable exceptions are as follows: 

A. In Radar 


24406 

Z$ESQ 

Earth's eccentricity squared 

24407 

Z1MESQ 

One minus earth's eccentricity squared 

24411 

Z$F 

Reciprocal of earth's radius - feet 

In Drag Table 


27763 

DRSC 

Negative earth radius - cm. Drag scalar 

Other Locations 


30546 

XRDT 

Velocity of vehicle with respect to 
reference body 

30535 

XRODT 

Unperturbed velocity of vehicle with 
respect to reference body 

30363 

T 

Time from launch - hrs 

30364 

h= T + 1 

Current integration interval - hrs 

303651 
30367 J 

\XI 

Perturbation displacement 

303701 
30372 J 

[■ XI+3 

Perturbation velocity 

303731 

30375-J 

j- D2XI 

Perturbation Acceleration 

50472 

TABLE 

Beginning of a block of locations described 
, in Section XI 
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MATHEMATICAL APPENDIX 


A. INTRODUCTION 

The problem of orbit determination over long time periods re- 
quires a precise technique for integrating the equations of motion. 
Reference 4 contains an analysis of integration procedure that yields 
the minimum loss of information due to the accumulation of numer- 
ical round-off errors. The Encke perturbation method has been 
shown to require minimum machine computation time for a mini- 
mum loss of numerical accuracy. The orbit prediction scheme 
presented herein uses a modified form of the Encke method with 
the initial position and velocity vectors replacing the conventional 
P and Q vectors of the Encke scheme. By avoiding reference to 
the position of perigee, it is possible to avoid numerical ambiguities 
arising from near circular orbits and orbits of low inclination. 
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B. EQUATIONS OF MOTION 


In a Newtonian system, the equations of motion of a particle in 
the gravitational field of n attracting bodies and subject to other 
perturbing accelerations such as thrust, drag, oblateness, radia- 
tion pressure, etc. are given by 



(B.l) 


These equations are put into observable form by referring them to 
a reference body c . The equations of motion of the reference body 
are 



Subtraction of Equation (B.2) from Equation (B.l) results in the 
equation of motion of the vehicle with respect to the reference 
body c. 


(B.2) 


R 


= ~ (/*v + Ac)tT ~ 


i=l 

ite 


R . 

VI 

R ci 

r 3 

VI 

r 3 

C 1 


L F. 


(B.3) 
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c . 


METHOD OF INTEGRATION 


If Equation (B. 3) is integrated directly by some numerical scheme, 
there results, after a number of step-by-step integrations, an accumu- 
lation of error which leads to inaccurate results. To avoid this loss in 
precision, it is convenient to write Equation (B. 3) in the form 

R vc = R^ + AR (C. 1) 

The velocity and displacement vectors can be written as 


R = 
vc 

R = 
vc 

The reference body (the one in whose sphere of influence the vehicle 
travels) is chosen so as to minimize the perturbations. 


Rk + AR 

R. + AR 
k 


(C.2) 
(C. 3) 


In this method 




is taken as 


and 


\ 


t = 


- (H + n ) 

V c 
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i^c 


R . 
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(C.4) 


(C.5) 


Equations (C. 4) constitute the equations of motion of the Kepler problem 
and are solved as described in Appendix D. 
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Equations (C. 5) are integrated numerically. The integration scheme - 
employed by the ITEM program is a sixth order backward difference 
scheme, initiated by a Runge-Kutta scheme. The routine used is 
"RW DE6F Floating Point Cowell Second Sum Integration of Second 
Order Differential Equations." (SHARE distribution #775.) 

As derived in the following subsection, the solution of the Kepler 

problem may be represented by the vectors R q , R q , the scalar a , 

and the rectification time t , 

o 

The rectification process consists of moving R vc , R vc into the 

locations R and R , t into t and computing a and n. For 
o o o 

computational convenience, the coefficients appearing in Equations (D.2) 
are also computed during rectification. 
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D. SOLUTION OF THE KEPLER TWO-BODY PROBLEM 


The unified formulation of the two-body problem is used for both 


elliptic and hyperbolic cases. 


J8 = 4\A • e 
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where 


d = R * R 
o o o 


r = ,8 F + r F + 




8 F. 


R = fR + g R 
o ° o 


R = f R +g R 
o o 




o 



a is determined from the modified Kepler equation 

7m At = 8 3 F i +r 0 8F 3 + “^ 2 F 2 

See Figure 4 for the two-body orbit which results from the solution of 
Equation (C. 4) with the initial conditions: 


R. (t ) = R (t ) = R 
k o ve' o' o 

k (t ) = R (t ) = R 
K O VC o o 


(D. 3) 


(D.4) 
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E. 


COMPUTATION OF PERTURBATION TERMS 


The terms accounting for the Encke term and the planetary per- 
turbations appearing on the right hand side of Equation (C. 5) involve 

R R o 

numerous terms of the form — — where R and R may differ 

o o O 

r r 

o 

only by small amounts. For the Encke term, for instance R - R^= £ 

which is small, and for the planetary perturbations, the difference is 

R which also often is small, 
ve 


A computation scheme, which avoids loss of precision due to the 

subtraction of nearly equal terms and which also is correct when R 

vc 

is not small, is employed. This scheme is described below: Find 



u = ~(R + ^ AR) • AR 
2 o 2 
r 
o 



o 


AR 

3 

r 

o 


R( u 3 + 3u 2 + 3u ) 


(-* 5 ) 


r 

o 


(E.l) 
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CONCLUSIONS 


The method presented yields accurate trajectories using relatively 
little computer time. Summarizing some of the important features: 

1. All significant solar system bodies may be included 
without undue complications. 

2. Since the perturbations only are integrated, the allow- 
able integration interval is fairly large over most of 
the path. Even in the vicinity of earth or another planet 
a relatively large interval (compared to other schemes) 
may be used without limiting the stability and accuracy 
of the solutions. 

3. The perturbations are kept small in two ways. First, 
the two-body orbit is rectified whenever the perturba- 
tions exceed a specified maximum value compared to 
the corresponding unperturbed values. This limits 
error build-up with respect to a particular reference 
body . Second, the reference body of the two-body 
problem is changed from earth, to sun, to planet ac- 
cordingly, as that reference body would contribute the 
largest perturbing force otherwise. 

4. This method will handle circular orbits, zero inclina- 
tion, etc. The problem is defined in terms of parameters 
which have real physical significance (namely, the posi- 
tion and velocity vectors) which are directly relatable to 
measurable quantities. 



G. OBLATENESS TERMS 


The oblateness perturbation terms in Equations (C. 5) are derived 
from the potential given by the following equation: 



(G. 1) 


where a^ is the equatorial radius of the earth. 


F = grad <p 


This vector can be written in the form 
F = l R + mk 


where k is a unit vector in the z-direction. 

a 2 


r 




(G.2) 
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+ ( r ) J 4oL 2 ( r* ' 


r35 ,z 3 15 z 

2 r 


<f)]} 


(0-3) 


The perturbation acceleration due to tesseral harmonics are computed. 


2 1 3 2 1 3 

Constants: J £ , J 3 , J g , X £ , X , X g 


J = -1.9E-6 

Li 


2 o 

X 2 = -21° 


J = -1.51E-6 

O 


1 

X 3 = ° 


(G.4) 


, 149E-6 


v 3 o 

X 3 = 22, 8 


In initialization: 


J 3 j 

C. = J. cos jX. 
-1 * ** * 


j 3 j 

S. = J. sin j X. 
1 1 ^ 1 


(G.5) 


x , y , z are x’ , y* , z* 


F 22 Term 


a 5 . 2 2, 

A 22 = ~ 2 X ' y > 


xy 

B 22=- 5 T 

r 


(G. 6) 
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„ 3fJ, f 2 2 1 

F 22x' = Tl C 2 (2x - XA 22 )+2S 2< y - xB 22 ) } 

F 22y' = ^5{ C 2 2( - 2y - yA 22 ) + 2S 2 2<X - yB 22 ) } 

„ 3fl f 2 2 ■, 

F 22 Z '= -1 C 2<- ZA 22 )+2S 2<- ZB 22>) 
r 

x, y , z are earth fixed coordinates, x, y in the equatorial 
plane, the x-axis toward the Greenwich meridian. Nutation is 
neglected. 

F Term 

OX 

^ ,a 2 2 2 , 

K 31 = (4z -x -y ) 

A = — K 
A 31 2 K 31 

r 

7y 

B 31=-i K 31 

r 

F 31x’ = I ^{ C 3 [K 31‘ 2x2-XA 31 ] + S 3 C " 2xy - XB 31 ] } 

F 31y‘ = I ^{ C 3 [ - 2xy - yA 31 ] + S 3 [K 31- 2y2 - yB 31 ] } 

F 31z' = I -7 { C 3 1[ SXZ - Z A 31 ] + S 3 C 8y * - Z B 31 ] } 
r 


(G.7) 


(G. 8) 


(G. 9) 
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F 33 Term 


7x 2 2 

A 33 = 1< x _3y > 
r 


7y 2 2 

B 33= ^< 3 * - y > 


{c 3 3 [ 3(x 2 - y 2 ) - XA 33 ] - S 3 [ 6xy-xB 33 ]} 

= ^{CgC-Sxy-yAgg] + S 3 3 [3(x 2 -y 2 ) -yB 33 ]} 

3, 



15 fi 

F 33x 

7 


r 


15 n 

F 33y’ 

7 

r 

"R 1 — 

15 jU 

F 33z' 

7 


{ C 3 '■" zA 33^ " S 3 ^ zB 33^1 


II 

F 22x* + F 31x’ + 

F 33x' 

n 

to* 

F 22y* + F 31y' + 

F 33y 

hj 

II 

F 22z’ + F 31 z’ + 

F 33z 


(G. 10) 


(G. 11) 


(G. 12) 
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H. TRANSFORMATION EQUATIONS FROM GEODETIC POLAR 
COORDINATES TO CARTESIAN COORDINATES* 


The geodetic polar coordinates in the program are referred to 
an ellipsoid of revolution. The equation of a cross section is given 
by 



(H.l) 


where 


b 2 - a 2 (l - e 2 ) 


The slope of the normal, along which h is measured is given by 


1 , a 2 z 

tan<£ •= — ; — = _ 1 " (See Figure 5) (H.2) 

dz h z x 

dx 


and 


tan <£>* 




t an 4> 


(l - e 2 ) tan 0 


Eliminating x between equations (H.l) and (H.2) and solving for z 
results in: 


a (l ~ e 2 ) sin <p 
(l - e 2 sin 2 4>) 1/2 


*For geocentric (i.e. e 2 “ 0) polar coordinates, c ~ s - 1. In this case the latitude 
input is interpreted as declination. 


H-l 



z 



Figure S. Relation Between Declination, Geocentric 
and Geodetic Latitudes 


and from equation (H.2) then 


a cos 0 
(l - e 2 sin 2 
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.In units of a , R and R are then given by equation (H.3) 

e 


c 

(l ~ e 2 sin 2 <fi) ™ 1/2 


s - 

(l-e 2 )c 


X = 

( c + h ) cos <fi cos ( 9 - < 9 Q ) 


y = 

( c + h) cos 0 sin (<9 ~& 0 ) 


z - 

( s + h ) s in <fi 


X = 

v £( s in 7 cos <fi - cos 7 cos A s in <56) 

cos (1 9 - 9 


- cos 7 sin A sin (<9 “ 0 O )| 


y - 

v £( sin 7 cos 4 > - cos 7 cos A sin <fi) 

s in (e - & 0 ) 


+ cos 7 sin A cos (d - 6 Q ) j 


z = 

v j~sin 7 sin <fc + cos 7 cos A cos 4 > | 



These equations include the effect of the rotation of the earth. The 
longitude of the vernal equinox (<9 0 ) at launch time is computed by 
the program from Newcomb's formula. 
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I. TRANSFORMATION EQUATIONS FOR RADAR SIMULATION 


The program computes sight angles (in an azimuth- elevation 
system), slant range and range rate data for up to 30 radar stations. 
The vehicle coordinates are transformed from a system of geocen- 
tric cartesian coordinates ( xyz ), the x-axis in the direction of the 
vernal equinox and the x-y plane in the equatorial plane of the earth 
to the required topocentric azimuth elevation system. This is ac- 
complished by a series of coordinate transformations as follows: 

1. A rotation of the coordinate system about the z-axis through 
an angle RA s so that x-y plane is in the meridian plane of the 
station. 


x ' - x cos RA + y sin RA 

S S 


y' - - x sin RA^ + y eos RA g 


i _ 

2 “ Z 


( 1 . 1 ) 


The velocity transformation must take the rotational velocity of the 
new coordinate system into account. 


x ' - y 1 co + x cos RA + y s in RA 

J e s' 7 s 


y ; - “x # oJ e x sinRA s + y cos RA g 

* * „ * 

z - z 


where x ' , y ' , z ' are the rotated coordinates and RA is the right 
ascension of the station and a> e is the sidereal rate of the earth's 
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rotation. The G.H.A. necessary to obtain RA s from the station 
longitude is computed by the program. 

2. A translation of the origin of the coordinate system from the 
center of the earth to the station in question 


x" - x' - (c + h)cos0 


tl 


y 


/ 

y 


z" - z 1 “ ( s + h ) s i ncfi 


( 1 . 2 ) 


where 


c - (l ~ e 2 sin 2 1/2 
s ~ (l - e 2 ) c 


• // 


y 


• « 


z 


z 


where x", y" , z" are the translated coordinates. <p is the geodetic 
latitude and h the height above sea level of the station in question. 

3. A rotation of (90 - <p) about the y" axis to place the ( x", z" ) 
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plane into the horizon plane 


x" # - x" sin0 + z" cos 0 


y 


at 


y 


II 


- “ x" cos 0 + z" sin 0 


9 at — 9 n ■ / , / 

x - x s m 0 + z cos 0 


( 1 . 3 ) 



z" = ~x" cos cp + z" sin<£ 


Now x'", y'", z"'are the coordinates of the vehicle in a topocentric 
azimuth elevation system, with z'" axis pointing to zenith and the 
x'" pointing south along the meridian. Range, range rate, azimuth 
and elevation are then given by 


p = (x'" 2 + y'" 2 + z'" 2 ) 1/2 = Slant range (1.4) 


x'" x'" + y'" y'" + z"' z'" 


P 


Range rate 


( 1 . 5 ) 


E - tan 


- 1 


(x'" 2 + y"' 2 ) 1/2 


= Elevation 


A' 


tan 


-ii 


X 


f 77 - A' A ' < 7t1 

(377 - A' A'>77j 


( 1 . 6 ) 


( 1 . 7 ) 
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J. TRIAXIAL MOON 


Triaxial lunar potential constants (as used in the ITEM Program) 

1 . The values of the constants A , B and C for the perturbation accelera- 
tions due to the triaxial moon may be calculated using data from the NASA 
earth model meeting. These constants are currently being used in the ITEM 
Program. 

2. The perturbation accelerations due to the triaxial moon are given by 
the partial derivatives of 






1 - 3 



(J.l) 


where 


u a 2 / 3 I_ \ 

m / C \ 

3 a 2 \2 ma 2 ) 

e ' A ■ m / 


A =' 



(J.2) 


B = 



J-l 



< 9 <£ 

<?x 


3xC 

—r-v 


d±_ _ _3yC 6CBy 
dy r s F r s 


(J.3) 


dcfi 3zC 6ACz 


F - 


dz r S *’ r S 


where 


F = h (—j- - 1) + B 


(¥ ■ ■)} 


Based on the NASA earth model meeting, the moments of inertia about the 
principal axes of the moon are*. 


I A = .88746 x 10 3s kg meters 2 


I B = .88764 x 10 35 kg meters 2 


I c = .88801 x 1 0 35 kg meters 2 

(See figure 6) 


Other constants are: 


/* = 19.9094165 x — -x~ 

nr 
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• Longest radius 


| (to earth) 

^-0 

\) 

Figure 6, Tri axial Moon 


(ER) 3 

r\ a A r\ f\ 

81.335 " • 

2447829 x 0 

nr 

(1738.09) 2 


(6378.165) 

- .0742595 (earth radii) 2 

5.975 x 10 27 

grams 




m 


m e _ 5.975 x 10 24 
81.335 " 81.335 

{1738. 09) 2 = 3.0209568 


= 7.34616 x 10 22 kg 

x 10 6 km 2 


The constants A , B and C may be calculated 


A 


B 


C 



( .88801 -.88746) 
.88801 


( .887 64 -.88 746) 
.88801 


619.36 x 10" 6 


202.70 x 10" 6 



For the ITEM Program the units of C are 

(earth radii) 5 


therefore 



, 2447829 ) { . 0742595 


3( .88801) x 10 29 


2 (7. 34616 x 10 22 ){ 3. 0209568 x 10 8 


C - 36.366998 x 10 4 earth mass x (earth radius) 2 
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In summary the constants used in the ITEM Program as based upon the 
NASA earth model meeting are: 


A = 619.36 x 1(T 6 


B = 202.70 x 10" 6 


C - 36.366998 x 10 4 earth mass x (earth radius) 2 
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K. DRAG COMPUTATION 


The drag acceleration is computed, assuming a spherically symmetric 
atmosphere rotating with the earth. Thus: 


where 


D = ~~2 | V e f f V e f f 


AR r 


JD 

m 


Veff ® OJ x R 


(K.l) 


co is the sidereal rotation rate vector of the earth. 
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L. COMPUTATION OF SUBSATELLITE POINT 


The geodetic coordinates of the subsatellite point are computed by the 
following method: 

The geocentric latitude (declination) is obtained from 


s m cp - — 


This latitude is then corrected to geodetic latitude by the formula 


<£ - + a 2 sin 2 cf>' + a. sin 4 4>' + a g sin 6 cp' + a 8 sin 8 4>‘ 


where 


a 2 " 1024 


§47 {512e 2 + I28e 4 + 60e 6 + 35e 8 } 


( e6 + e8 } "^ r { 4e6 + 3e8 ) 


a 4 = ■ TOM? {64e 4 + 48e 6 + 35e 8 } 


1 r. 4 15e 8 e 8 

_ , \4e 4 + 2e 6 + e 8 / + — - — 

K- 2 >. J oec.3 tc. 


256r 3 16r< 


^ (4e 6 + 5e 8 } - ^ {e 8 + e 8 } 


{4e 6 + 3e 8 } 
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e 8 f 5 64 252 320 1 

a 8 “ 2048 Y r + r 2 “ r 3 + 7 T J 

e = the eccentricity of the earth 


r = the distance from earth's center 


See Reference 5. 

The geodetic height is then given by 


h = r cos - ]/l - e 2 sin 2 


The longitude is obtained by subtracting the sidereal time of 
Greenwich from the right ascension given by 


tan RA = 


(L.4) 


(L.5) 
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M. POLAR COORDINATES REFERRED TO THE MOON 


Moon longitude and latitude are defined in the coordinate system 
described in appendix J. If the vehicle coordinates with respect to 
the center of the moon in this coordinate system are given by x , y, 
z , r , then 


6 - longitude = tan 1 


4 > - latitude 


sm 


-l. 
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N. SHADOW LOGIC 


A coordinate system is set up in the plane defined by the centers 
of the light-emitting source, the shadowing body, and the probe. Both 
bodies are assumed to be spherical, and hence all testing can be car- 
ried out in this plane. The diagram in Figure 7 shows this plane. 

The coordinates are defined by unit vectors i and j: 

R ot . . . . 

J_ = “ — ; J_'_l = 1 5 JL’J = 0 

cl 

lj_ = -di+R vc ; * * 

where 

d - R • i 
— vc — 


Vehicle coordinates in this system are given by: 


x = 
v 



d 


y = R • j 
J v —vc — 


= c. 


,2 x 2 -i 

■ d + r j 
vc 


1/2 


z = R • K = 0 
v — vc — 


(N. 1) 
(N.2) 

(N.3) 

(N.4) 

(N.5) 
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•1. Shadow Parameters 


a) The tips of the umbra and penumbra cones are: 





d = 
P 



b) The slopes of the bounding lines are: 


r 

c 


sin a - cos 0 = _ 
u u d 


u 


*- cos a = sm 0 - ; 1 -{ — ^ 
u u L \d / - 
u 


sin a 


tana = 


u 


u cos a 


u 


sin 0 


tan 0 = 


u 


u cos 0 


u 


cos 0 = sin a = — 
P P d 


r 2 - 

■ Q r 1 ' C ^ 12 

sin 0 = cos a =i 1 - s — 1 ’ 
p p L \d / j 

P 

sin a 


tana = — 1 ■— c 
p cos a 


sin 0 


tan 0 = 


p cos 0 


c) Refraction Correction: (UMBRA) 

a' = a - € , 0' = 0 - € 

u u u u 


sin a' 
u 


tan a' 
u 


tan 0' 
u 


sin a cos € - cos a sin € 
u u 

tan a - tan € 
u 


1 + tan a tan € 
u 

tan 0 - tan € 
u 


1 + tan 0 tan € 
u 
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d' = 
u 


s in a 


u 


d) Refraction Correction: (PENUMBRA) 


Both 

e < a 

e> a 

; a = a 


P 

P 

P P 

sin a' 

= 1 sin a 

COS € - 

cos a sin e 1 

P 

P 


P 


tana 

- tan e 



tan a' = 
P 


tan 0' = 
P 


1 + tan a tan e 
P 

tan 0 - tan e 
P 

1 + tan 0 tan e 
P 


d’ = - sign (tana' ) ■ : 

p & ' p ’ sin a' 


The equations of the bounding lines are given below. 


2. The Testing Procedure 


0 Line 
P 


% = 


tan 0' 
P 


v 

tan 0' 
P 


x s 0 
v 


x < 0 

V 


Sunlight 


Go to next test 
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y 1 

- ( x - d' ) tan a.' > 

0 

Sunlight 

V 

v P p 


y v' 

- ( x - d’ ) tan a' = 
v p p 

0 

Sunlight penumbra boundary 


< 

0 

Go to next test 



s 0 Penumbra 
< 0 Go to next test 


> 0 Penumbra 

= 1y | - ( x^ - ) tan = 0 Shadow penumbra boundary 

< 0 Shadow 


and are stored and saved. The crossing times 
are found by linearly interpolating for 0-values of Q and 
respectively, to guarantee that crossing from one region 
into another always occurs across these boundaries. 
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O. SOLAR RADIATION PRESSURE 


The radiation pressure subroutine computes the force of solar radiation 
on the spacecraft if an appropriate pressure coefficient is used. The calcu- 
lation relies on the shadow routine to set a trigger to multiply the pressure 
coefficient by 1.0, 0.5, or 0.0 for full sunlight, penumbra or umbra, respec- 
tively. Therefore, the shadow subroutine must be used in conjunction with 
the radiation pressure routine for most cases. If the spacecraft is known 
to be continually in sunlight the number 1, may be loaded into SHADIN and 
thus elaborate shadow testing may be avoided. 


C p. AR vs 


RP 


m r 


vs 


( 0 . 1 ) 


(See Section VIII-H-4. for definition of symbols.) 

This radiation pressure subroutine has been found to be inexact for 
satellites of large area to mass ratio since it only controls the pressure to 
the nearest integration step. For such spacecraft (e.g., balloons) several 
degrees error in true anomaly may result after 100 days unless the integra- 
tion is carried exactly to the boundaries. A modification to achieve this 
increased precision is available and will be included in future versions of 
the program. 
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P. ECLIPTIC COORDINATES 


The ecliptic coordinates are an approximate set obtained by a simple 
rotation of the equatorial coordinates about the x-axis through a fixed angle 
i = 23°26’31" which is approximately the true obliquity for Jan 0.0, 1962. 
More exact coordinates may be obtained by changing NE , (unit normal to 
the ecliptic) as desired. 
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Q. MOON ROTATING AND FIXED COORDINATE SYSTEM 


Geocentric coordinates of the vehicle based on the earth-moon 
plane are generated from the geocentric equatorial radius vector 
to the vehicle, R yE , the geocentric unit vector in the direction of 
the moon R me , and the vector in the direction of the moon's velocity, 

^me • 

Coordinates in the rotating system, XROT etc., are found by using 
the current values of these vectors at each time step in the relations. 


XROT 


R 


VE 



“or = «ve -(«*e* *«) 


ZROT 


R 


VE 


H 


ME 


(Q.i) 


where 



^me x ^ME 

i^ME x 


For the fixed axis system XINJ, etc., the initial vectors R ME (t 0 ) 
and R me (t 0 ) at the time of injection are used with the current 
value of R ve . 
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TRAJECTORY SEARCH 


The program provides a search routine to obtain selected trajec- 
tories. The search is based on linear theory and varies the polar load 
input quantities (independent variables) to search for desired dependent 
variables. There are twelve possible dependent variables to select 
from, although a maximum of seven of the twelve may be used in any 
given search. The quantities, at present, are 


1, O, oo, t (pericenter time), and 
P 


Tp (pericenter radius) 


The above five variables at the moon and at the earth (in the 
case of earth return trajectories) constitute the first ten variables. 
(These quantities are normally referred to the equatorial plane. How- 
ever, the earth-moon plane is also available. See Section VIH.) 


In addition, the components of the impact parameter vector 
( B * T , B • R ) may be selected. They are referred to the ecliptic 
plane for Mars and Venus trajectories and the moon's orbital plane 
for lunar trajectories . The number of independent variables need 
not be equal to the number of dependent variables for this routine 
to operate . (See Section VHI-C-1 . ) 

This search routine is time consuming if the initial conditions 
are poorly approximated. Before using this routine, two things should 
be done. 
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1. A first guess of the initial conditions of the nominal trajec- 
tory should be obtained from a patched conic or a similar 
search program. 

2. The number of variables should be kept to a minimum. It 
is planned to automate the iteration scheme to go from two- 
body, to patched conic, to full trajectory, and to increase 
the number of variables to be adjusted, in optimal fashion. 
Even in its present form, however, it is extremely useful. 


The iterator uses a modified version of the MIN-MAX Principle 
(Reference 6). 

A.. is the matrix of partial s 

Ax. is the vector of changes in the independent 
variables 

X„ is a diagonal matrix of weights 
y is a vector of residuals 

The system to be solved is 


(A.. A.. + X.. )Ax. = A., y. 
' ji i] n ' i Ji i 


A.. + X... = 
i] 


B.. 

13 


A.. 

Ji 


z. 

l 


e> 


R-2 



Procedure 


The system 

B.. Ax. = z. 
i] i 1 

is solved. If the value of Ax. is greater than SIZER , some arbi- 
trary amount, set 

B.. = B..+ X.. 

13 ij ii 

and solve the system again. Repeat these operations until Ax. is 
less than or equal to SIZER. Now, run a new nominal trajectory 
with the new independent variable 

x. = x. + Ax. 
i i a 

a) If the new residuals y are greater than the previous ones, set 

B.. = B.. + X.. 
i] i] 11 

Solve the system again and continue solving until the new 
residuals are less than the old. Now the system is ready 
for a new iteration. 

b) If the new residuals y. are less than the old, set 

B.. = B.. - X.. 
ij i] ii 

Solve the system again and continue solving until the new 
residuals are greater or equal to the old. 

The iteration continues until either the maximum number of 
iterations (input) is exceeded or the residuals are less than or equal 
to an input tolerance. 
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S. EQUATIO NS FOR FLIGHT PATH AZIM UT H AND FLIGHT PATH ANGLE 


A subroutine computes the flight path azimuth and flight path angle 
with the following equations: 

1. Flight path angle 


7 = 



A 

•N 


(S.l) 


N is the vertical unit vector. In the geodetic system $ is given by 


N = 


cos (p cos (< 9 - 5 0 ) , cos ^ sin [d - & 0 ) , sin<p 


In the geocentric system 4> is replaced by <£', Alternatively, in the latter 
system 


N 



2. Flight path azimuth 


A sin 


- 1 


" 1 

fy 

cos y 

tv 


cos 




sin 




A 


cos 


-1 


1 fi . 

cos 7 cos <£ j"v" “ s in 7 s m 0 


(S.2) 


Both formulas are used to determine the proper quadrant of A . To obtain 
the geocentric output, e 2 = 0, <p is replaced by declination S .= <p‘ . 
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T. OSCULATING ELEMENTS 


The osculating elements are obtained from the following equations: 


a 



(T.l) 


n = fj} /2 \a\~ 3/2 


(T.2) 


e cos E 
e cosh E 



(T.3) 


e 

e 


sin E 
sinh E 



(T.4) 


f E - e s in E 

lesinhE-E (T.5) 


t 


p 



(T.6) 


The angles fi , co , i are obtained from the vectors H and P, where 


H = R x R 


(T.7) 


eP 



R - 



(T .8) 


In terms of these vectors: 


H Z 

cos i = “j^ - in the first or fourth quadrant (T.9) 
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sin ft 


COS ft - 


COS Cp 


s in co 


H 

X 

h sin i 



h sin i 


= P cos ft + P sin ft 

x y 


p 

z 

sini 


(T.10) 


(T .11) 
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U. IMPACT PARAMETERS 


The "impact parameters" are coordinates in the "impact" plane. This 
plane passes thru the body (T-planet or the moon) and is normal to the 
incoming asymptote. The direction cosines of the asymptote are given by 
equations (U.l, U.2) in terms of unit vectors P (Appendix T) and 

Q 


s 

A A A 

In the plane defined by S as its normal, two unit vectors T IMP and R IMp are 
defined. T IMP is parallel to the ecliptic plane for mars and venus impacts 
and to the moon's orbital plane for moon impacts. Explicitly 


H 6 

h x F 


= T (p - /F-T)q) 


(U.l) 

(U.2) 


T 


IMP 


N X S 

In x si 


(U.3) 


where N is the unit normal to the ecliptic plane, or the moon's orbital 
plane. R IMp is normal to both § and T IMp . B IMp is the vector from the body 
to the vehicle as it crosses the impact plane. The information computed 
are the dot products 


®imp ‘ ^imp anc * ®imp ‘ ^ 


IMP 
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W. EQUATIONS FOR TRANS LUNAR PLANE INPUT 


The translunar plane input is designed to permit easy visualiza- 
tion of the geometric relationships between initial conditions for 
circumlunar trajectories and the motion of the moon. (See Figure 1). 

The initial conditions are given in a coordinate system referred 
to the translunar plane. This system has its x-axis along the ascend- 
ing node of the vehicle with respect to the moon' s orbital plane, its 
y-axis in the translunar plane, at right angles to the ascending node, 
in the direction of motion. In this coordinate system, initial position 
and velocity vectors are given by 


(r B + h)cos ¥ 
(r B + h)sin¥ 
0 


Here r B is the radius of the body of departure (earth or moon). 

x tl ~ v sin (y - ¥ ) 


y tl - v cos ( 7 ~ ’P ) 


The translunar plane is positioned by giving its inclination i TL 
with respect to the moon' s orbital plane and the lunar lead angle <p, 
the angle between the moon's position at injection and the descend- 
ing node. The vectors R TL and R TL may then be transformed into 
the equatorial system by the following series of rotations: 


X TL 

y TL = 

Z TL 


(W.l) 


(W.2) 
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1. A rotation - i TL about the x TL axis will rotate the translunar 
plane into the moon's orbital plane. 

2. A rotation of tt ~ (\. M + 4 ) about the new z-axis will refer 
the moon's orbital plane coordinate system to the ascending node 
of the moon's orbital plane (with respect to the equator) as x-axis. 

Here A_ M stands for the argument of latitude of the moon. These 
rotations are performed by multiplying R TL and R TL by the matrix: 


-cos (Ajj +4) 

sin +< p) 

- sin + 4 ) sin i TL 

-sin (>•« +4) 

- cos ( \j + 4) cos i TL 

cos {^+4) sin i TL 

0 

sin i TL 

cos i TL 


(W-3) 


3. The moon's orbital plane (MOP) is rotated about its node 
through an angle - i M (the inclination of the MOP). 

4. The ascending node is brought into coincidence with the 
vernal equinox by a rotation - 0 M , These two rotations are embodied 
in the matrix 


/ cos fi M - s in n M cos i M s in 0 M s in i M 
! - sinO M cosO M cosi M -cos ^ M sin i M 

\ 0 sin i M cos 

and thus: 


(W-4) 


R — ( BA ) R tl 


R = ( BA ) R tl 


(W.5) 
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Y. CHANGE OF INDEPENDENT VARIABLE - BETA MODE 


According to the standard Encke method, we introduce a differential 
equation 


P 



(Y.l) 


In the construction of the closed-form solution for (Y. 1), a parameter 
0 arises. It is related to t by Kepler's equation, 

f(0) 

t = t + -p=J- (Y. 2) 

° fj. 

where f is a transcendental function of /3 and is obtained by summing 
several power series. 

If t is taken as the independent variable, Equation (Y. 2) has to 
be solved for 0 by an iterative method, requiring numerous time- 
consuming evaluations of the function f for each integration step. 

Using /3 as the independent variable, however, only requires a single 
evaluation. 


It remains, of course, to see what becomes of Equation (Y. 1) and 


l 


= x - p 

= -p(r 



+ F 


(Y. 3) 
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if 0 is the independent variable. We have, from Kepler's equation, 
that 


dt_ _ IpI 


at any point along the solution of (Y. 1). Thus 

,</£" , . IpI 

p * p rrr and p - p -7= 

1 p 1 V£T 


at any point along the solution of (Y. 1) and the initial conditions become 

x | x | 

P<0 O ) = x o and p'(3 o ) = 

when 


0 = fl(t ) = 0 
o o 


Now the solution for (Y. 1), p and pV, can be written in closed form 
for any 0. As auxiliary quantities in this solution, we have jp | and 

p • p' 

D = — — - . They are computed as functions of 0 before p and p' are 

IPI 

known; that is, with accuracy at least as good as that of p and p\ 

Not only are they needed and easy to compute, but they also have the 
interesting property that 

dt_ _ IpI 

and 

d 2 t = _D_ 

d0 2 7]T 

Thus Equation (Y. 1) is solved more economically in terms of 0 than 
in terms of t. 


(Y.4) 


(Y.5) 
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- Now we turn to Equation (Y. 3). To treat it, we want to express 
£" in terms of \ . From (Y. 5) we have that 


f i t dt f 1 Pi 
5 " 


(Y. 6) 


Differentiating with respect to /3 , 


t" ^ \ — 1 . t ' 


'• VJ- 




, D 

P I 'Tp" 


= € 


fiL + *v-£, 
* IpI 


= - Ip 


2 / x 


f( 


Y i |3 
x 


,3/ M 


F + e’ri 


(Y.7) 


Thus (Y. 7) is the equation to be integrated numerically, instead of 

I P P D 

(Y. 3). The coefficients and : — - can be calculated with much 

IpI... ... 

more accuracy than the factors involving £ , since they depend only 
on the two-body solution. For analysis of error propagation, we 
write (Y.7) as 



(Y. 8) 
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The mechanics of the procedure, then, are easy to enumerate. The 

initial conditions are x and x . Let 

o o 


P(t ) = x 
o o 


p'(t o ) = 


X X 

o o 


Using these initial conditions, evaluate t 
value of 0 to be considered. 


1 PP D. 

m ’ \py 


p , p for each 


(Y.9) 


Let 4 q = 0. Using these initial conditions, integrate Equation 

(Y. 7) to get 4(0) and 4'(0). Note that the first two terms on the 
right-hand side of Equation (Y. 7) are functions of x and possibly x’ . 
These are obtained by 


x<0) - p(, 8) +4(0) 

x'(/3) = p'(/3) +4’(0) 


(Y. 10) 


If at any point x is required, it can be found from 


x[t(0)] = x’(0) 


I p<0) I 


(Y. 11) 


Depending on the rectification control logic, there will be places 

where the solution to Equation (Y. 1) must be started over. At this 

point, the values t , x, x become the new t , x , and x , while 
c o o o 

0 , 4 » and 4’ a.re reset to zero. 
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Comparison of Modes 

a) It is immediately apparent that eliminating the necessity 
of iteratively solving Equation (Y. 2) will substantially 
increase the speed of computation. 

b) An important advantage arises farther from eliminating 
the sometimes ponderous logic which supplies initial 
guesses for the iterative process and guarantees con- 
vergence of the solution. 

c) A third advantage of the 0 -method is not quite so apparent, 
but no less important. It is well known that the size of the 
integration time step can be increased as the distance from 
the center of attraction increases. This change of the inte- 
gration interval requires a, cumbersome restart procedure. 
An examination of Equation (Y. 5) shows that equal intervals 
of correspond to time intervals of increasing lengths as 
the distance increases. The time interval thus automatically 
expands and contracts correctly without outside intervention. 

d) Geometric stopping and printing conditions can usually be 
conveniently expressed in terms of j3 , whereas they often 
require iterative determinations of the time. This advan- 
tage, however, is slight. 


Y-5 



e) If state vectors are required at fixed times, an iteration is 
necessary to find the corresponding value of jS. In this case 
the j3 -method is no better, and no worse, than the standard 
methods. If such vectors are required at frequent, closely 
spaced, time points (as in orbit determination, for instance), 
the advantage of the j3 -method is marginal. 
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